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Energy transfer through a dissociated diatomic gas in 
Couette flow 


By JOHN F. CLARKE 
The College of Aeronautics, Cranfield 


(Received 2 April 1958) 


SUMMARY 
The transfer of energy through a dissociated diatomic gas in 

Couette flow is considered, taking oxygen as a numerical example. 

The two extremes of chemical equilibrium flow and chemically 

frozen flow are dealt with in detail, and it is shown that the surface 

reaction rate is of prime importance in the latter case. The 

chemical rate equations in the gas phase are used to estimate the 

probable chemical state of the gas mixture, this being deduced 

from the ratio of a characteristic chemical reaction time to a 

characteristic time for atom diffusion across the layer. The 

influence of the surface reaction appears to spread outwards 

through the flow from the wall as gas-phase chemical reaction 

times decrease. For practical values of the surface reaction rate 

on a metallic wall, the energy transfer rate may be significantly 

lower in chemically frozen flow than in chemical equilibrium 

flow under otherwise similar circumstances. 

Similar phenomena to those discussed will arise in the more 

complicated case of boundary layer flows, so that a treatment of 

the simpler type of shear layer represented by Couette flow may be 

of some value in assessing the relative importance of the various 

parameters. 

1. INTRODUCTION 

There is at present much interest in the flow of gases at temperatures 
much greater than those which can be withstood by all known solids and, 
consequently, it is of some importance to know the rate at which heat will 
be transferred from the gas to the walls of its container, or those of some 
body passing through it. In heating the gas to such high temperatures 
it is highly probable that its polyatomic components will dissociate and 
absorb a large amount of energy in the process. The dissociation energy 
then exists as a potential source of heat energy, which may be released 
following a suitable chemical reaction. The usual heat transfer problem 
is therefore complicated by the need to consider the rates at which such 
chemical reactions may occur. In this connection it is necessary to consider 
both homogeneous reactions (those which occur in the gas phase) and 
heterogeneous reactions, namely those which occur at a solid surface 
adjacent to the gas. 
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The general problem is therefore one of some complexity although 
practical solutions have been found in some cases. ‘The most comprehensive 
of these to date is the work of Fay & Riddell (1958), who confine their 
attention to the forward stagnation region of blunt-nosed bodies. Fay & 
Riddell’s solution represents a more exact analysis of the boundary-layer 
equations, following along the lines of Lees’ (1956) largely heuristic 
discussion of the same problem. Both references take explicit account of 
the difference in the rate of transfer of heat energy by thermal conduction 
and the rate of transfer of dissociation energy, the transport medium in 
the latter case being simple mass diffusion of atoms. It has been assumed 
in both cases that air behaves as a single diatomic gas. Previous solutions 
of the dissociation problem have been given by Moore (1952), Crown (1952) 
and Beckwith (1953). The first reference deals with the flat plate problem 
whilst the second two deal with forward stagnation regions, but all assume 
that the gas is in chemical equilibrium and imply that thermal energy and 
mass diffuse at identical rates (more strictly, they assume equality of 
the appropriate diffusion coefficients). 

The complexity of the problem is increased in all of these works by the 
need to consider the boundary layer on a flat plate or at a stagnation point, 
giving rise to the familiar non-linear ditferential equations of boundary 
layer theory. However, much that is characteristic of the boundary layer 
is reproduced in the simpler viscous flow known as Couette flow. The 
lack of variation of properties in the direction of flow then creates notable 
simplifications, a fact which has been exploited by Illingworth (1950) for 
the case of a fluid of constant chemical composition. Under this 
restriction, Illingworth’s work represents the only solution of the full 
Navier-Stokes equations valid for all Revnolds and Mach numbers. More 
recently Liepmann & Bleviss (1956) have used Couette flow to exhibit 
some of the properties of a shear layer of dissociated gas. ‘Their treatment 
implies equality of mass and thermal energy diffusion coefficients, however, 
and only considers flow at chemical equilibrium, but it extends the gas 
temperatures into the range for which ionization becomes appreciable. 

In the following sections we consider the question of energy transfer 
through a simple dissociated diatomic gas in Couette flow (using some values 
appropriate to oxygen for numerical illustrations). ‘The energy equation 
is readily integrated, yielding a relation between enthalpy, atom concentra- 
tion, and velocity. A simplification of the treatment is achieved following 
the assumption that the gas behaves as an ideal dissociating gas, a concept 
introduced by Lighthill (1956). Real gases such as oxygen and nitrogen 
appear to be very nearly ideal in this sense in a temperature range sufficient 
to cover problems ot practical interest and a simple expression for enthalpy 
in terms of temperature and concentration results. From here on the 
problem resolves itself into the determination of atom concentration in 
terms of either temperature or velocity in any particular set of circumstances. 
Once this is achieved, the skin friction coefficient and Stanton, or Nusselt, 
number follows directly from integration of the momentum equations. 
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In order to determine the atom concentration it is necessary to consider 
the kinetics of the dissociation and recombination reactions. While not 
difficult in principle, the treatment results in an unwieldy non-linear 
differential equation which would require numerical integration. Results 
are given below for two extreme cases, namely very fast and very slow 
homogeneous reaction rates, the full rate equations being used only to give 
an estimate of the conditions under which such assumptions may plausibly 
be made. The heterogeneous, surface reaction is considered in both cases 
and is shown to have a controlling influence on energy transfer rate when 
gas-phase reactions occur very slowly. 

It may be permissible to interpolate here some brief remarks on the 
mechanisms involved in the energy transfer process. All three of the 
transport phenomena are explicitly concerned in the process, the basic 
mechanism in the gas phase being thermal conduction. ‘Thermal energy 
is the sum of the energies contained in the translational, rotational, and 
vibrational modes of the gas molecules, so that its rate of transport must 
be a function of the facility with which these energies can be interchanged. 
When such interchange occurs so rapidly that it is reasonable to assume that 
each internal mode is in equilibrium at the local temperature, the classical 
Eucken correction is derived which modifies the monatomic thermal 
conductivity value to account for the extra degrees of freedom of polyatomic 
molecules. At the temperatures which we shall consider here, the Eucken 
correction is probably accurate enough, but, in any case, it will be found 
that thermal conductivity can be absorbed into certain dimensionless 
numbers which may reasonably be assumed constant across the gas laver. 
Energy transfer also arises from the dissipation of ordered kinetic energy 
of mass motion into random thermal energy due to viscous actions in the 
gas. The final rate of heat recovery by these means is affected by the ratio 
of the vorticity and heat energy diffusion coefficients, in other words, the 
Prandtl number. 

When the gas consists of a mixture of molecules and their dissociated 
atoms a further potential source of heat energy exists. Each atom in the 
mixture can be visualized as a vehicle for the transport of dissociation energy 
as it diffuses through the gas layer, the energy being released in the form 
of heat as a result of a recombination reaction. The rate of heat recovery 
in this case will be found to be a function of the ratio of the mass diffusion 
coefficient and the thermal energy diffusion coefficient (the Lewis number) 
in exact analogy with the viscous dissipation phenomenon. Furthermore, 
this heat recovery rate depends on the appropriate chemical reaction rates 
and it will also depend on whether the dissociation energy is released in 
the gas phase or directly at the solid surface. 


2. ‘THE CONSERVATION EQUATIONS 
The mass conservation equation for the 7th species of a general mixture 
of reacting gases can be expressed in the form 
div[7,(q+q,;)] = ,, (1) 
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where n; is the number density of species 7 and q; is its diffusion velocity 
vector. q is the mass average velocity vector and w,; is the ‘source’ term 
which results from production of species 7 due to chemical reactions, the 
units of w; being number of species 7 produced per unit volume per unit 
time. If m; is the mass of a particle of the 7th species, the general mass 
conservation requirement is }m;w;= 0. Since Smj,n;q; = 0 as a result 
of the definition of diffusion velocity, multiplication of equation (1) by m; 
followed by summation over all 7 leads to 
div pq = V, (2) 

where p = }mjn,, the mixture density. Equation (2) is the usual form 
of the continuity relation. 

The momentum balance is unchanged when chemical reactions occur 
and is written as 


p(q.grad)q+divp = 0, (3) 


where p is the pressure tensor (see, e.g. Hirschfelder et al. 1954, eqns. 
(7.2-23) et seq. for definition). 

The energy equation can be derived by considering a particular 
elemental mass of gas moving with the mass average velocity. In a steady 
flow this gains internal energy EF at a rate pq.grad FE units of energy per 
unit volume per unit time due to convection. £ is the internal energy 
of the gas mixture per unit mass and is equal to p-!}m;n;e,, where e; is 
the internal energy of a particle of the 7th species above some arbitrarily 
determined zero level. For example, for a binary mixture of symmetrical 
diatomic molecules and their dissociated atoms we can take the internal 
energy of a molecule, e,,, as zero at the absolute zero of temperature. Then 
the internal energy of each atom, e,, at the same temperature will be 
one-half of the energy required to dissociate the molecule at this temperature. 
The elemental mass under consideration will also gain kinetic energy of 
mass motion by convection, but this is accounted for by the action of the 
stresses exerted by the fluid on the element. In addition these stresses 
do an amount of work —pdivq in compressing the element (p is thermo- 
dynamic pressure) and dissipate an amount of mechanical energy ® which 
appears in the form of heat. @ is the familiar dissipation function. 
Finally, the element gains energy from two other sources, namely from 
thermal conduction and by diffusion of the separate species across its 
boundaries. ‘Thermal conduction provides a term div(Agrad 7), where T 
is temperature and A is the coefficient of thermal conductivity, whilst 
ditfusion results in a rate of gain of energy per unit volume per unit time 
equal to —div(}n;/;q;). hf; is the enthalpy of a particle of species 7 and 
appears here in preference to internal energy e; since each component gas 
in the mixture must perform flow work on the gas of its own species which 
is in front of it. It should be noted that the additional kinetic energy of 
each species which arises from the difference between mass average and 
ditfusion velocities has been neglected here, but this is an assumption 
which is generally made (e.g. Chapman & Cowling 1939, p. 145). The 
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contributions from thermal conduction and mass diffusion can be combined 
in a single quantity qy, the energy flux vector, defined as 
qo = —Agrad T+ Yn; h,q,, (4) 
so that the energy conservation requirement can ultimately be written in 
the form 
pq. grad E = —divq,—pdivq+®. (5) 
We shall neglect the thermal diffusion contribution to diffusion velocity 
since this is generally negligible, particularly when we have to deal with 
large concentration gradients. (A description of thermal diffusion is 
given by von Karman (1956).) In that case the diffusion velocity of species i 
is given by 
2 
q; = = dm, D,; grad(n;/n) (6) 
pn; j 
(see, e.g. Hirschfelder et al. 1954, eqns. (7.4-3)). The quantity 7 is the total 
particle density while ; and n; are the particle densities of the 7th and 
jth species. m,; is the mass of a particle of the jth species. Dj; is the 
diffusion coefficient for diffusion of species 7 through species j, the 
summation being taken over all species except 7. We shall confine attention 
to binary gas mixtures, in which case (6) reduces to 


9 
n? ‘ 
Gq, = — M,Z, grad(n,/n), (7) 
pny 
where %,, is now a binary diffusion coefficient. Some simplification of (7) 
results from the definition of a mass ratio, c,, such that 
pc; = m;N;. (8) 
Plainly }c; = 1 and (7) then reduces to 
9; = —LZ grade, (9) 
with a similar equation for qs. (Note Zs, = Zo.) 

The conservation equations can also be written in terms of mass 
fractions, as can the energy flux vector qy. We may also use enthalpy 
per unit mass rather than per particle, but, before making these modifications 
we will simplify the equations to apply to Couette flow. 


3. COUETTE FLOW 

The variation of all properties in the x-direction is zero in Couette 
flow (see figure 1 for definitions). In order to satisfy the continuity 
equation, equation (2), pv = const., « being the mass average velocity 
normal to the walls. But « = 0 at y = 0 and at y = 6; whence it must be 
zero everywhere. 

With these simplifications and treating only a binary mixture consisting 
of symmetrical diatomic molecules with their dissociated atoms, the conser- 
vation equations are as follows. Equation (1) gives 


£ (pea) = 64 M (10) 
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for continuity of the atomic species. (We now write suffix a to indicate 
atom, sufhx m for molecule.) 7, is the atom diffusion velocity normal to 
the walls. The momentum equations (3) give 


df{ du dp 
— —_ — 0 =—, 
7 (u ) dy it) 
and the energy equation (5) becomes 
d{.dT d du\* 
—( A— ) — — (pc, Hy qt plm Am Um) + ul =) = 9. 12 
z( zB) dy (pc, a a pe m mt a o() ( ) 


H, is enthalpy of the ith species per unit mass and yp is the dynamic 
coefficient of viscosity. Making use of equation (9) for the diffusion 
velocities and the fact that c,+c,, = 1, equation (12) can be rewritten 





=f ar’ d dc, es 2 . 
Kasil (0 Vitter eee | Dottie aw b «= O, Ki 
. du zy) 7 PZam dy e =) C) 
LY 
Upper wall | Suffix §” 
— with | / 
veocityU~__ || / UU 
oe 
i.U § 
| Xx 
Stationary / Suffix ‘w’ 
Lower wall. 


Figure 1. Couette flow. 


where the ditference H,—H,, has been set equal to D, the dissociation 
energy per unit mass at the absolute zero of temperature. This is an 
approximation which implies equality of the specific heats (per unit mass) 
of atomic and molecular species. In the present circumstances this is 
quite a reasonable approximation and will be dealt with further in § + below. 
Equation (11) can be integrated at once to give p = const. and 
i = 7 = Const, = +, (14) 
dy 

(7 is the shear stress, suffix w implying evaluation at the lower wall, y = 0). 
Putting equation (14) into equation (13) permits integration of the latter 
relation, giving 


1T Ic 
pple oD am D==4 + ur7,,, = const. (1 
dy dy 


wm 
— 
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Evaluating the constant in equation (15) when y = 0, it is readily seen 
with the aid of equation (+) to be equal to —q,,, the energy flux into the 
lower wall. The total enthalpy of the mixture, H, is given by 


H=c,H,+c, fi, 
so that, since p is constant here, 
-« @f dH dc 
ee ay ey a pe ¢ 
" dy dy dy v ‘. 
where 
C.=c,C,,+e 


p a pa m pm 


(16a) 


C,,. and C,,,, are the specific heats of the atomic and molecular species at 
constant pressure, so that C,, is the specific heat of the mixture when its 
chemical composition is frozen. ‘The approximation D ~ H,—H,, has 
been made in equation (16). Then equation (15) can be written as 


\ dH , A \dc . du? a 
= —- + D(p A tn — ee ial oe e —_ = — Gy (17) 
C, dy C,/ dy 2 dy 
Defining Prandtl number Pr and Lewis number Le as 
Pr = uC, A, le = pC, Gia, (18) 
respectively, and assuming that they are constant, integration of equation (17) 
between y = 0 and y = y gives 
H—-H,,+ D(Le—1)(c,—¢,,.) + $Pru? = —uPrq,,/7,,- (19) 
This result follows from the fact that 
4 dy “ 1 dy u 
;— = | -—_— du — ws 
lt pe du T 


u being zero at y= 0. When y=6, then H =H,, c,=c,, andu=U; 
equation (19) therefore gives 


H, — H,,+ (Le-—1)D(c,, —¢,,,.) + $Pr U? = —U Pr q,,/7,,- (20) 


we 


The recovery enthalpy, H,, is defined as the value of H,, when q,, = 0. 
Whence 


H, = H,+ D(Le—-1)(c,, —¢,,) + $Pr U?, (21) 
where c,, 1s the concentration of atoms at y = 0 under zero heat transfer 
conditions. Equation (20) can now be rewritten 


—(q,,/7,,)U Pr = H,-H,,+ D(Le-1)(¢,,-Cyw)s (22) 
so that, defining a skin friction coefficient, C,, and a Stanton number, St, 
such that 

ty, = dp, UFC,, (23) 


— 4, = py U\H, —H,, + D(Le—1)(Car—Caw)} Sts (24) 


equation (22) shows that 
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This is the form of Reynolds analogy for Couette flow of a dissociating gas 
mixture, the analogy factor being exactly 1/Pr in this case. This result 
is also given by Liepmann & Bleviss (1956) but is shown to be perfectly 
general here (provided St and C;, are defined as in equations (23) and (24)), 
since no mention has so far been made about the chemical state of the gas. 

The problem is now reduced to that of finding 7, (or C,). Integration 
of equation (14) gives 


U 


Wy, = | BP du, 
or, in dimensionless form, 
pu’ m 
y’Re,.C, =2| —du', (26) 
; JG ty 


where y’ = y/6, u’ = u/U and Re,, = p,, Ud/p,,. 

The viscosity . is dependent upon temperature and chemical composition, 
the value of C, deduced from equation (26) depending rather critically 
on this variation. There is some evidence, however, to show that p for 
the type of gas mixtures which we are considering continues to obey 
Sutherlands law quite well up to very high temperatures (e.g. 9000° K), 
and with this in mind we shall assume p» proportional to 71*. While this 
assumption may not be strictly justifiable, it simplifies the analysis and 
should enable reasonable comparisons to be made under varying conditions. 
Noting that Nusselt number, Nu = StPrRe,., equation (25) shows 
that we now can express equation (26) as 

-u “1 
vy Nu= | VT’ du’; inparticular Nu = | \/T’ du, (27) 

' “0 

where 7” = T/T,,. Equations (19) and (22) show that 

H' —H/ + D'(Le-1)(c, —¢,,) + byPrM?2 u? 

= u'\H — H+ D'(Le—1)(Cup—Caw)}, (28) 
where H’, H’, H’ and D’ are H, H,,, H, and D divided by RT, IV. 
R is the universal gas constant and W’,, is the molecular weight of the 
molecules. .W? = U?W,,/(yRT,,.) and is similar to a Mach number, based 
on the speed of sound at the lower plate. This is only strictly true if the 
gas is undissociated at the lower wall, y then being the ratio of specific heats, 
but it is best regarded as simply a convenient non-dimensional group In 
the present circumstances. ‘The enthalpy is a function of temperature and 
concentration which is known a priori, so that, once c, is established as 
a function of either 7” or wu’ (or both), equations (27) and (28) enable us 
to find Nw and hence velocity, temperature and concentration profiles. 

Evaluation of concentration now constitutes the major problem and 
to this end one is led to consider the kinetics of the dissociation and 


recombination reactions. Before proceeding with a discussion of this 
topic we will introduce the concept of the ideal dissociating gas, which 
is due to Lighthill (1956), and which permits some further simplifications 
in the treatment of the problem. 
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4. ‘THE IDEAL DISSOCIATING GAS 
It can be shown that the number ratio of the constituents of a dissociating 
diatomic gas in equilibrium is given by 


tac — Ja exp(—d/kT) (29) 
—, - 
where sufhx e implies equilibrium value. The quantity d is the dissociation 
energy per molecule and k is Boltzmann’s constant. f, and f,, are partition 
functions for the atomic and molecular species, and may be taken to be 
products of the simple partition functions of the appropriate degrees of 
freedom, provided cross-couplings between different freedoms are negligible 
(see Moelwyn-Hughes 1947). The mixture density is p = m,n,+2m,n,, 
and it is readily shown that equation (29) gives 
a — Plexp(—d/kT), (30) 
l-—Cue =p 
where p, = m, f2/2f,,: this quantity is a characteristic density (Lighthill 
1956), which is found to be roughly constant for oxygen and nitrogen 
over the temperature range 1000°K to 7000° K, being about 150 gm/cc 
for O, and 130 gm/cc for Nj. The thermal equation of state for the mixture 
is 
p = (1+c¢,)p(R/W,, )T, (31) 
assuming that each component gas in the mixture is thermally perfect, 
so that equation (30) could be rewritten in the form 
ce. 6—-spg RT 
1-c2, pW, 
This is a rather more convenient form for present purposes. 

The approximation p, = const. requires f? to be proportional to f,, 
and it can be shown (Lighthill 1956) that this is equivalent to fixing the 
amount of energy stored in the vibrational and rotational modes of the 
molecules at 3k7T. (In fact this energy will vary from kT to 2kT, depending 
on the degree to which the vibrational mode is excited.) ‘The enthalpies 


exp(— W,, D/RT). (32) 





per unit mass of each species are then 

H_,, = 5RT/W,,+D, H,, = 4RT/W,.> (33) 
and so H,—H,, = D+ RT/W,,. However, RT/W,, is much smaller than D 
for the temperatures which are of interest here and H,,—H,, ~ D is a very 
reasonable approximation. (In fact, for a real gas, oxygen for example, 
the vibrational mode of the molecules is very nearly fully excited at such 
is nearer to 4(RT/W,,), so that the approximation 


nm 


temperatures and H,, 
is better than would be suggested by treating oxygen as an ideal dissociating 
gas.) 
The specific enthalpy of the whole mixture is plainly 
H = (4+c,)RT/W,, +¢, D. (34) 
Equation (34) is used later to eliminate the H’ quantities from equation (28), 
(note that H’ = (4+c,)T’+c,D’). If it is found reasonable to assume 
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that the mixture is close to a chemical equilibrium state (i.e. ¢, ~ ¢,,), 
equation (32) can be used to eliminate c, from equation (28) in terms of 
temperature. Otherwise, c, must be found from the chemical rate 
equations. 


5. REACTION KINETICS 
Consider the dissociation reaction of the symmetrical diatomic species, 
A,. The stoichiometric equation describing the reaction is 


Ry : = 
A,+X@2A+X (35) 
k 


, 
in the usual notation (e.g. Penner 1955). The specific reaction rate constants 
for forward and reverse reactions, k, and k,, are generally expressed in 
terms of the concentrations of the respective reactants measured in 
moles/cc, for reactions of order higher than the first. We shall write (A,) 
to mean ‘concentration of species A, in moles/cc’ and similarly for the 
other species. ‘The species X in equation (35) may be either A, or A in 
a pure mixture, so that (X) represents the total number of moles of mixture 
per cc. Accordingly (X) = p/RT, where p is thermodynamic pressure 
and R the universal gas constant, so long as the individual gases A, and A 
can be treated as thermally perfect. 

The net rate of production of the atomic species A, in moles/cc per sec, 


1s 
HAY 9 - 7 2 
alee 2R,(A2)(X) — 2k (A)?(X), 
or, in terms of mole fractions, 
( (A) r\o. 7\3 2 
. a 2k,(X)*x,, — 2k,(X)*x,. (36) 


{The mole fractions are x (A,)/(X), x, = (A) (X).) Under chemical 
equilibrium conditions, d(A)/dt = 0 and equation (36) gives 


m 


(37) 


K., is the equilibrium constant in terms of mole fractions (or what is 
equivalent, number density ratios) and can be found from equation (29). 
Thus, knowing either k, or k,, the other can be found from equation (37). 
Let us now write equation (36) in the form 
d , > \< 9 Ax, 
— [x,(X)] = —2k,(X)8[2x,, — x2, + Ax, — x,,Ax,]—, 
dt saa 
aey and assume Av, <1. The reaction will now proceed 
at very nearly constant p and 7, so that (X) and x,,, may be assumed constant 
to all intents and purposes. Then 
1 dAx,, Fe 
Ax. <a = 


where Ax, = x,-% 


2x =e 
| ~ const., (38) 


1—x,, 





24S) * 
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and we can define a characteristic time, ¢,, such that 


1 1 dAx, 


rs _ Ax, dt 
Equation (38) shows that ¢, is constant under the stated conditions; it 
is in fact the time taken for a small deviation of atom concentration to fall 
to 1/e of its equilibrium value at conditions p and 7. Thus t, can be used 
as a measure of the rate at which the chemical reaction proceeds in the 
direction of equilibrium, and can be compared with some characteristic 
flow time in order to give an indication of the probable chemical state of the 
gas. The obvious choice of flow parameter is the time taken for an atom 
to diffuse across the layer. In Couette flow equation (9) gives 


(39) 


Gi ik, 
d ¢, dy’ 





and we can reasonably assume (c7! dc,/dy’) to be of order unity. The 
binary diffusion coefficient, Y,,,,, 1s of order 1 cm?/sec at room temperatures 
(~ 300° K) and 1 atmosphere pressure and varies roughly as T*?/p. Thus 
v,, ~ (T/300)??(p(atm.)5)—! so that a characteristic time for diffusion of an 
atom across the layer, ¢,, may be defined as 


Pa TF \ 3/2 
Dg tte (550 ss (40) 
t, re) 300 po- 


( 





Equations (38), (39) and (40) now give 
t, a, 1-x,, (= "ft T \32 1 
t, 2|2x,,— Xz, p / R,\300/  p(atm.)d?’ 


fe E bas |e vi (p in atmospheres). (41) 


2X 4¢— X2,_| p87, 





or 





We can say that if ¢,/t, < 1 then the flow will effectively be in chemical 
equilibrium in the gas phase, whilst if ¢,/t, > 1 then chemical reactions 
occur so slowly in the gas phase that chemical composition will, effectively, 
not vary on this account. Plainly the factor in square brackets in 
equation (41) can exert considerable influence on the ratio of the 
characteristic times through its influence on ¢,. The mole fraction x,, 
is never either zero or unity in a finite temperature range (see equation (29)), 
although it may approach these values extremely closely. Equations (38) 
and (39) merely indicate that any excess or deficiency of atoms above or 
below the relevant equilibrium value is adjusted extremely rapidly when the 
mixture is predominantly of the atomic species, whilst the reverse is true 
when the gas is predominantly composed of molecules. The reasons are 
that recombination is a three-body process (two atoms plus any other particle) 
which is therefore relatively improbable when few atoms exist. Dissociation 
is a two-body process requiring, however, a certain high energy level 
encounter between the particles concerned and this is a relatively improbable 
occurrence at the lower temperatures where the mixture is mainly composed 
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of molecules. Within the range 0:1 < x,, < 0-9 the quantity in question 
varies by less than 10? in magnitude. 

Before equation (41) can be used to provide an indication of the 
magnitude of ¢,/t, one must know the magnitude and temperature 
dependence of k,. Here the available evidence is rather indecisive. 
Feldman (1957) suggests that a lower limit for k, in air (treated as a 
single diatomic gas) is about 10!* cc?/mole?/sec for temperatures up to 
about 5000° K, the trend being for k, to decrease with rising temperature. 
Fay & Riddell (1958) write the reverse specific reaction rate constant as 
k, = K, T-*?, where K, is a constant, and estimate K, to be 5 x 10'(300)?? 
for oxygen; in other words, 


k, =~ 5 x 10'4(7T/300)-3? cc?/mole?/sec. 
Hirschfelder (1956), however, quotes a value of 
k, =~ 10'6(7/300)>? cc?/mole?/sec 


for oxygen, basing this result on the theoretical work of Eyring, Gershinowitz 
& Sun (1935) for recombination of atomic hydrogen, whilst Camm & Keck 
(1957), by monitoring the radiative relaxation zone behind a moving 
normal shock wave, deduce a value k, ~ 10!'® cc?/mole?/sec in pure oxygen 
or nitrogen, which is substantially invariant with temperature. In view of 
this conflicting evidence it does not seem unreasonable for the present to 
set k, = 10" cc?/mole?/sec and ignore its temperature variation. Apart 
from Fay & Riddell’s, the estimates quoted above lie within roughly an 
order of magnitude of this value either way in the interesting temperature 
range for dissociation. ‘Therefore we shall write 


a ts ~ 10-5 ps? (42) 


and, remembering that this estimate may vary by one or two orders of 
magnitude in either direction, use equation (42) to estimate the probable 
chemical state of the gas mixture. 

[n particular it can be seen that the flow must be very close to a chemical 
equilibrium state for pressures of 1 atmosphere or greater, so long as 6 is 
not much less than 1 cm, but that a decrease of pressure, say down to 
10-*-10-* atmospheres, combined with shear layer thicknesses around 
1 mm will rapidly lead towards the freezing of homogeneous reactions*. 

So far we have not considered the influence of surface reactions on the 
flow pattern. When the ¢,¢, ratio is very small it seems reasonable to 
suppose that particles remain in a particular flow region for a sufficient 
time to allow the chemical equilibrium state to be very closely approached 
there. Consequently, this state would constitute a good first approximation 
to the true state of affairs and any iterative procedure founded upon it 
would be expected to converge very rapidly. Under these conditions, the 


* In applying considerations of this sort to an object flying at hypersonic speeds, 
it must be remembered that pressures in the boundary layer on the object will usually 
be many times the undisturbed, atmospheric value. 
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species’ concentrations and mass fluxes become functions of temperature 
and pressure which will not admit the imposition of arbitrary boundary 
values at the walls y = 0 and y = 6. The surface reaction rate which would 
specify these boundary conditions will vary very widely (for different wall 
materials, for example) so that it may seem that conditions of chemical 
equilibrium in the gas phase and arbitrarily fast surface reactions are 
incompatible. 

For the case of heat transfer in a stagnant gas, Hirschfelder (1956) 
has shown that the two conditions are reconcilable following the introduction 
of transients in the flow variables which have appreciable values only in 
a thin gas layer adjacent to the surface. The methods used by Hirschfelder, 
to evaluate thethickness of this layer are valid for Couette flow. The example 
quoted there, for energy transfer through oxygen at one atmosphere pressure 
with no surface reaction, gives the thickness of the layer as less than 10-3 cm, 
the actual thickness decreasing with increase of wall temperature (N.B. 
temperatures lie in the range 2000-6000° K). For low wall temperature 
values, the fast homogeneous reaction rates cause the atom concentration 
to fall to zero (for all practical purposes) some distance from the wall and 
the problem of surface reactions does not arise. It is interesting to note 
from equation (42) that with p= 1 atmosphere and 6 < 10% cm, ¢,/tg 
becomes greater than unity, a fact which would lead one to suspect the 
equilibrium approximation. Hirschfelder shows that for plate separation 
distances greater than the distortion layer thickness, energy transfer rate 
is virtually unaffected by heterogeneous reaction rate. 

The continuity requirement for the atomic species can be expressed 
in the following form, using equations (10) and (36) and writing the 
diffusion velocity v, in terms of concentration gradient, 


a 


d de ; ; p \8 4c? , fl—-c - 
=P gc =o Ph he = a ae Sh I. 4: 
Leva |= & Walzer) Lasep &(TGe) | 


Kis the equilibrium constant in terms of mole fractions and can be found 
from equations (29) and (37), 








n= 4RT 
—_ ae i ae >K = - T 7 44 
MN ne pu a —_ ? ) ( ) 


where 7', = d/k, the characteristic dissociation temperature. For oxygen 
T,=59000° K and p,=150gm/cc; whence pK, = 1540T exp(— 59 000/7), 
with p measured in atmospheres. Writing k, as K,(T)’ cc?/mole?/sec, 
equation (43) for oxygen becomes 


d dc a 4¢2 
ot oP =f f = 29» 10 °K, OT | ee = 
alee | orl oe 


~ 1540 exp(— 5900/7) ~(7—*2) | (45) 
p \lt+e, 
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p being measured in atmospheres. One could now eliminate 7 from 
equation (45) in terms of c, and u and the y-derivatives could be replaced 
by u-derivatives, using equations (28) and (27). This would result in a 
differential equation for c, in terms of u, but involving the Nusselt number, 
which is not known a priori. ‘The solution would not be seriously influenced 
by the value of Nw and no doubt some plausible value for this quantity 
could be employed at the start of an iteration scheme. However, the 
equation is very unwieldy, and no attempt has been made to solve it here. 
Instead we shall go on to consider the two extreme cases of t,/t, < 1 and 
t./t,, > 1, assuming in the first case that atom concentration is equal to 
the local equilibrium value, and in the second case that no gas phase reactions 
occur at all. In the latter case it becomes necessary to consider the details 


of the surface reaction. 


6. CHEMICAL EQUILIBRIUM FLOW 

When ¢,/t,, < 1 it can be assumed that c, ~ c¢,, everywhere at the local 
temperature and pressure. It should be observed that this does not imply 
that the net rate of production of atoms is, or even approaches, zero. Indeed 
the reverse is true; the rate of production of a particular species can be 
considered to be ‘very fast’ in these circumstances. It is just that the 
chemical reaction can proceed at the required rate with concentrations 
not too far removed from an equilibrium value. Mathematically, it could 
be said that the coefficient of the term in square brackets on the right-hand 
side of equation (43) approaches « while c, — c,, in such a way as to cause 
the product to approach a suitable finite, non-zero limit. Then equation (43) 
is not required and c,, can be found from equation (32). Putting this value 
into equation (28) provides the required (7", 1’) relation. 

Since we do not seek great accuracy it is simplest to solve equation (28) 
graphically and, likewise, to integrate equation (27) graphically to find Nu. 
As a numerical example, some values have been obtained in this way for 
energy transfer through oxygen at one atmosphere pressure. The 
temperature at the upper plate has been taken to be 5000 K, the lower 
plate being maintained at 7',.= 1000 K. The dimensionless dissociation 
energy per unit mass, D’, is therefore equal to 59 in this case, and 
Py Py» = 0°39 x 108, The Prandtl number has been taken to be 0-75, while 
y at the lower wall is equal to 4 for the ideal dissociating gas. These values 
are reasonable approximations 1n the circumstances as we are more interested 
in trends and comparisons than in actual values. The (7",w’) relation 
has been found under these conditions for two ./,. values, namely ./,,. = 2 
and W,.= 20, and several Lewis number values between 1-0 and 2:0. 
The value of Le for oxygen is probably about 1-4, the other values being 
included to demonstrate the effects of mass diffusion on the profiles, etc. 

Dimensionless velocity, temperature, concentration of atoms and 
enthalpy profiles are shown in figures 2 to 5 for MW, = 2 and figures 6 to 9 
for W,.= 20. Figure 10 shows the variation of Nusselt number with 


Lewis number for the two .W/,,. values. 
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All profiles are, to a greater or lesser extent, affected by the value of the 
Lewis number, although the variations are so small in the case of velocity 
distribution when 1/,,. = 2 that they do not show up on the scale of figure 2. 
Similarly Nu depends slightly on Le, although the variation is not 
detectable when W,.= 2. The results of §3 show that the rate of energy 
transfer into the lower wall is given by 
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Figure 2. Velocity profile in chemical equi- Figure 3. Temperature profiles in chemical 
librium flow. M,, = 2, T,, = 1000° K, equilibrium flow M,, = 2, T,, = 1000° K, 
T; 5000° K, p 1 atmosphere. All T, = 5000 K, p = 1 atmosphere. 
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Figure 4. Atom concentration profiles in Figure 5. Enthalpy profiles in chemical 
chemical equilibrium flow. ©,, = 2, equilibrium flow. M,, = 2, T,, = 1000° K, 
T,. = 1000°K, 7; = 5000°K, p= 1 T; = 5000° K, p = 1 atmosphere. 
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Figure 6. Velocity profiles in chemical Figure 7. Atom concentration profiles in 


equilibrium flow. MW, = 20, T,, = 1000° K, chemical equilibrium flow. M,, = 20, 
T; = 5000° K, p = 1 atmosphere. T,,, = 1000° K, T; = 5000°K,p = latmos- 
phere. 
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Figure 8. Temperature profiles in chemical equilibrium flow. 


M,, = 20, T,, = 1000° K, Ts = 5000° K, p = 


1 atmosphere. 
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Since Vu is substantially invariant with Le, the energy transfer rate 
is almost a linear function of (Le—1), other things being equal. For the 
examples evaluated here the ratio q,,(Le = 1-4)/q,.(Le = 1:0) is about 1-3 
when M,,.=2 and 1:1 when M,,=20. An analysis which assumes 
Le=1-0 may therefore give energy transfer rates which are quite 
significantly low, particularly for the case in which viscous dissipation 
contributes a relatively small term. (A practical example of such a flow 
occurs in the region of the forward stagnation point of hypersonic vehicles. ) 















































It should be observed that for all cases considered in this section, 
de, dy = de,, dy + O0as y > 0. Therefore all the energy is transferred into 
the lower wall by thermal conduction alone and temperature gradients 
there should be larger for larger values of Le. ‘This effect can be seen in 
figure 3, when V/,. = 2, but does not show up on figure 8. Inspection 
of the solution in the latter case shows that in fact the (7", v’) curves cross 
over near to 7’ = | and the picture is then similar to the lower half of 
figure 3. The effect of Lewis number on the profiles may be qualitatively 
explained as follows, viewing an increase in Lewis number as an increase 
in the ability of atoms to diffuse through the layer. Then for larger values 
of Le these particles can penetrate to regions closer to the wall before 
recombining and giving uv their dissociation energy as heat, resulting in 
the steeper temperature gradients there. in the outer parts of the layer, 
viscous heating tends to maintain the atom concentration at its upper wall 
value, indeed when .V/,. = 20 the viscous heating is to intense as to cause 
an increase in atom concentration (see figure 7; note c,, = 0-965 in this 
case) and rapidly leads to virtually complete dissociation. ‘The increased 
mobility of the atoms (implied by larger Le values), permits them to escape 
from this region more readily, however, and the viscous heating can only 
maintain a smaller zone of high concentration when Le is large. The net 
result is that a greater proportion of the viscous heating 1s available in the 
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form of thermal energy when mass ditfuses less readily, as exemplified by 
the very large temperature peak in figure 8 when Le = 1-0. The temperature 
values there will be modified somewhat by ionization of the oxygen atoms, 
a fact which may well cause Le to be of some importance in magneto- 
gasdynamic applications. 

The effects of variation of wall temperature ratio, 7, 7, on Nu are 
rather complicated, as can be appreciated from an inspection of equations 
(27) and (28). The rapid rise of ¢,, with temperature is an important 
factor and it has not been found possible to relate Vu to T, 7,,.in any simple 
manner. However, in an attempt to obtain a reasonable correlation between 
these quantities for the single case, 7. = 1000" K, it was found that better 
results were obtained if the peak value of temperature in the layer, Ras 
was used in place of 7,. This choice atfected the high ./,. values most 
since viscous heating gave 7), values in excess of 7, in all cases. Only for 
values of 7, approaching 7). was the effect felt at /7,.= 2. Some results 
of plotting Vu against C = (px), (pu), are shown in figure 11, suffix p 
indicating evaluation at peak temperature conditions. The general quantity 
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Figure 11. Variation of Nusselt number with wall temperature ratio, Ts; 7T,,, in 
chemical equilibrium flow 
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C is important in boundary-layer theory, where it appears as an explicit 
factor in the boundary-layer equations. Fay & Riddell (1958) have 
demonstrated its significance in the stagnation region where a factor 
[(ex); (pu),.]°* appears in the heat transfer rate value, sufhx 6 implying 
evaluation at the edge of the boundary laver in this case. It is not too 
surprising to tind that C correlates wall temperature variation quite well 
in Couette flow, therefore, although it must be stressed that the examples 
quoted here are for one lower wall temperature value only (7. = 1000° K). 
Since flow in a stagnation region is generally one for which viscous 


dissipation plays only a small part in determining energy transfer rate, 
it tends to be closer to Couette flow at low .V/,. values, where 7, 1s less 
than, or at most of the same order as, 7,. On this basis, wall temperature 
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ratio effects in boundary-layer flows with large viscous heating terms may 
perhaps be better correlated using a peak temperature in the layer, as is 
so for the W,, = 20 case in Couette flow. The small amount of evidence 
which has been derived here to support this suggestion is in no way 
conclusive, but may perhaps be borne in mind in future studies. 

The values of C used here are equal to [(1+c,,),(7,/T,,.)'7]“, since 
Caew IS zero for all practical purposes when 7. = 1000°K. It appears 
from figure 11 that Nw varies as some power of C which differs in different 
ranges. When 7, = 20, the lowest value of 7, is about 3800° K, so that 
a fair degree of dissociation still exists there (c,, ~ 0-47 at T = 3800° K, 
p = 1 atmosphere). The same is true of the first three values of C when 
M,, = 2 (see figure 11), but at C ~ 0-5 (where 7, ~ 3250° K and ¢,,, = 0-1) 
the slope alters abruptly. It must be concluded that peak temperature 
and degree of dissociation have a significant effect on Nusselt number. 


7. CHEMICALLY FROZEN FLOW 

Having considered some aspects of chemical equilibrium flow (t,/t,, < 1) 
attention will now be directed towards flows at the other extreme, for which 
t/t, > 1. In this case, chemical reactions occur so slowly in the gas phase 
that they can reasonably be neglected altogether in a first approximation 
and the term w, in equation (10) put equal to zero. This state will be called 
chemically frozen flow and §5 has indicated that it may be a valid concept 
at low pressures and/or small values of 6. 

To solve this problem it is necessary to consider the details of the surface 
reaction which may take place. ‘The heterogeneous recombination reaction 
occurs between a gas utom and an atom adsorbed (or chemisorbed) by the 
wall. Observations (fur example Shuler & Laidler (1949)) have shown this 
to be a first-order process, requiring the surface to be fully covered with 
adsorbed particles. The efficiency of the wall as a catalyst for the 
recombination reaction can be defined as the ratio of the number of gas 
atoms being converted into molecules to the total number of gas atoms 
striking the wall per unit area per unit time, and will be denoted by I. 
For the type of flow represented here by Couette flow, the microscopical 
system is not too far removed from its equilibrium state and for present 
purposes it can be assumed that the actual distribution of molecular 
velocities is Maxwellian. The Maxwellian distribution function, F, is 
expressed in terms of the peculiar velocities of the particles of a particular 
species (vector q’ with components uw’, v7; and w,), this velocity being defined 
as the velocity of random, thermal motion of a particle relative to the local 
mass average velocity. The mean of q; over all particles of species 7 is the 
diffusion velocity, q;. The distribution function F is given by 


\32 
. m ee - 
F=n ‘ee exp(—m,q*/2kT), (47) 
ZakT) “*P | 
where n; is the number density of species 7, m; is the mass of a particle and 


k is Boltzmann’s constant. The number of particles whose velocities lie 
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in the range q’ to q, + dq; is F du’ dv, dw;. Since the mass average velocity 
at the wall is zero in both Couette and boundary-layer flows, the number 
of particles of the th kind striking the wall per unit area per unit time will 


be 
px “x “0 - : 
Sim = | du, | dw; | UF dv, =n, RT (27m; kT,,)—*?, 
- -¢ . x 
or 
Sip = Xap (22m; kT,,)-™2, (48) 


where x,,, is the mole fraction (= 7;,./n,.) and p is the thermodynamic 
pressure (=17,,kT,,). Equation (48) now provides an estimate of the 
rate at which gas atoms strike the wall and from the definition of recombina- 
tion coefficient I’, the mass flux of atoms into the wall is seen to be I'm, s,,,.. 
With no reactions occurring in the gas phase, the continuity requirement 
(equation (10)) becomes simply 








o@, Ha E. ( Re ) pm, a (49) 
" dy l+c,,,.J/ (27m, kT,,)'? 
in particular, for oxygen, 
J dc, 350 ' C.- <0) 
is Sa (3 =) nie 


We will only deal in detail with the lower wall, assuming that conditions 
at y = 6 correspond to the equilibrium state at the appropriate temperature 
and pressure. ‘lhe quantity p%,,,, in equation (49) can be written as pu Sc, 
where Sc is the Schmidt number (Sc = Pr/Le). Sc is constant through 
the layer since both Le and Pr have been assumed constant, so that writing 
proportional to ,/7 as before, equation (49) assumes the dimensionless 


de, 1 dScUp 26m (51) 
dy’ T’ p,(27m,kT,,)!? \1 7a; 


(dc, dy’) (dc, du’ )(du’ dy’) = (dc, du')(.Nu \ T’) 


from equation (27), so that equation (31) can be written 


form 





But 


] 


where 4 is a constant. Integration of equation (52) yields 








“eae iC; -c,,,.)u' 4 Gres (53) 
since c, = c,, When wu’ 1 and c, =c,,. When u’ =0. The constant A is 
given by 

bS8cC ] 2c, . os 
, oy ere Bs! en end (54) 
Hy(2em, RT)? Nu\1+ cay 


The concentration of atoms at the lower wall can now be determined 
as a function of recombination coefficient in any particular circumstances, 
once the Nusselt number has been found. By writing the dimensionless 
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enthalpy values appearing in equation (28) in terms of 7” and c¢,, (see § 4), 
equation (53) can be used to eliminate c, in terms of vu’. Then the (7",w’) 
relation required for evaluation of Nu is found to be 
[(4+¢,5)(Ts/T,,.) — (4+ Caw) + byPrMi lu’ — dy PrMz uu”? 
(4+.6, 4) + (Cas — Coy )U’ 


T’,u’) relation is a function of ¢,,,. it would be possible to use 


T’ 





Since the ( 

equation (55) to find Vu for some chosen ¢,,,. value and then to use equation 
(54) to find the value of I’ appropriate to the chosen conditions. (Values 
of T,, T,,, M,,, p, ete. would be known since they are required to specify 
the problem.) However, it is found that the (7",u’) relation does not 
depend very critically on the value of c,,,, so that Nw is substantially 
invariant with this quantity for any chosen 7°,/T',, ¢,;, M,,, etc. Lewis 
number does not appear in equation (55), except through its influence 
on ¢,,,, Whence it can be inferred that Nu is independent of Le for all 
practical purposes. 

The energy flux into the lower wall is 


dT’ dc - 
<i ee +p,.G,.. Di =a 5¢ 
Lu 4a), 2 Pu Lam (3), . ( )) 


and is in part due to thermal conduction, in part due to atom diffusion 
with release of dissociation energy directly at the wall. Since the (7",u’) 
relation is independent of c,,,. it can be seen that the (7’,y’) relation is 
likewise independent of this quantity; therefore, for any given values of 
T;, T,,, 5, ¢,, and M,,, heat transfer by thermal conduction takes place 
at the same rate for all values of [in chemically frozen flow. In other words, 
the temperature profile through the layer does not depend on the surface 
reaction rate. The remaining mechanism of energy transfer, however, 
depends critically on the value of . The second term on the right-hand 
side of (56) can be rewritten as (p,, Z,,,, DNu/5)(dc,,/du’) since 7’ = 1 
when y’ = 0, and (53) shows that (dc,/du’),.- 9 = Cas—Caws Where ¢,,, is 
found as a function of I from (54). It can be seen from (54) that c,,,, is a 
minimum when I’ has its maximum value of unity and furthermore that 
¢,,, can never be equal to zero. (It is, of course, always less than c,, unless 
I’ is negative. ‘This would imply that the surface reaction was a dissociation 
rather than a recombination reaction, however, requiring 7',, values higher 
than would be practically possible.) From this result, the conclusion can 
be drawn that it is not possible to recover all the dissociation energy in 
chemically frozen flow that may otherwise be recovered in chemical 
equilibrium flow under similar conditions, provided that the 7’, value 
is not too high. (It seems unlikely that 7,,. could be allowed to go much 
above 1000-1500 K in a practical case.) 

For comparison with the previous, chemical equilibrium case some 
values have been calculated for 7, = 5000 K, T,,, = 1000° K and M,,. = 2, 
the gas being oxygen as before. ‘The value of c,, has been set equal to 1-0 
and the Schmidt number Sc = 0-75/1-4. A reasonable estimate of pw for 
the mixture is wz = 5 x 10-4(7,1000)'? gm/cmsec (being roughly the value 
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for molecular oxygen at high temperatures), bearing in mind the assumptions 
previously made with regard to this quantity. The Nusselt number, 
calculated graphically, has the value 1-76 and equation (54) gives 


10? _.. ens 
1 Ce = 730 (5 as), 


(In the event that [ = 0, equation (54) gives c,,,,. = « 
\z 





‘ag = 1-0, and Nu can 
be evaluated analytically using ( For the example quoted here, 
Nu = 1-73, in good agreement with the graphical evaluation for c,,,. < ¢,5-) 
Figures 12 and 13 compare the velocity and temperature profiles in the 
two cases (N.B. Le 
independent of 6 so long as it is not too small in the chemical equilibrium 


(5 
55 


1-4 for the equilibrium flow), the profiles being 


case. ‘The curves marked ‘chemically frozen flow’ are good approximations 
for all surface reaction The filling out of the equilibrium state 
temperature profile due to gas phase recombination reactions is apparent in 
figure 13. Figure 14 comnares the (‘ dimensionless ’) concentration profiles. 

When the product pd = 10-3, equation (57) gives 0-15(1—c2,,) = Tc 
and this relation is shown plotted in figure 15. 


rates. 


aw 


The least value of c in 


a 
this case is almost equal to 0-15. It can be seen that pd = 10~% is consistent 
with the requirement for chemically frozen flow, since equation (42) could 
be written as ¢, t,, ~ 10 p, and frozen flow would in all probability occur 
for p < 107! (In passing, it appears to be quite 
possible that frozen flow conditions may be realized on a model in a con- 


to 10-* atmospheres. 


ventional shock tube, provided the low pressure end can be sucked down 
sufficiently far. Such conditions may also arise during satellite re-entry 


into the earth’s atmosphere. ) 
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Figure 12. Comparison of velocity profiles. 


AM m2, 1000° K, T; = 5000° K. 
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Figure 13. Comparison of temperature profiles. 


M,, = 2,T, = 


1 


1000° K, 7; = 


5000° K. 


- 1-0, the rate of energy transfer is a maximum and it is 


interesting to compare it with the chemical equilibrium value in otherwise 
Denoting final energy transfer rate in chemically 
frozen flow by qg,,, and in chemical equilibrium flow by q,,,, the ratio 


similar circumstances. 
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Figure 14. Comparison of atom concentra- Figure 15. Surface recombination coefficient 
tion profiles. M,, = 2, T,, = 1000°K, vs atom concentration at the wall. 
T, = 5000 K. M,,=2, T, = 1000? K; T, = 5000° K. 
pd = 19-° atm. cm. 


Vuy/Iwe = 9-85 for the examples quoted here. Thus it would appear that 
the ‘fully catalytic wall’ is a less effective means in these particular circum- 
stances of recovering dissociation energy than gas phase recombination. 

Couette flow can be interpreted as a crude approximation to a laminar 
boundary-layer flow in which the effects of convection have been neglected 
compared with those of diffusion, the Couette flow thickness 5 being regarded 
as a local boundary-layer thickness. It is therefore interesting to compare 
the present theory with that of Fay & Riddell (1958), bearing in mind that 
actual numerical values may be widely different, whilst the general effects 
of the various parameters on the final result are in all probability similar. 

For the heat transfer rate, — q,,, in the stagnation region of a blunt-nosed 
body Fay & Riddell give 


where Re is a Reynolds number, Re = 





—9, = —— Nu,(H, + }U?—-H,,), 
al ; : ' 
ium. where x = distance measured round the nose from the stagnation point, 
suthx w refers to the wall, and sufhx 6 to the edge of the boundary layer. 
U is velocity at the edge of the boundary layer, so that H,+3U? = H,, 
the stagnation enthalpy. .Vu, is a Nusselt number found to i given by 
: o oh Delta ve” 
Nu, = 0-67Re!? ; (1+ A(Le—1)(D/HA, cus}; 
files. Pw he 


Uxp,,/u,, and De,, is the amount 


of energy per unit mass in the flow outside the boundary layer which is 


stored as dissociation energy. 


A is a constant almost equal to 0-5 which 


varies somewhat according as the flow is in chemical equilibrium or is 


chemically frozen. 


(In their later paper, 4 = 1 and the Lewis number 
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is raised to a power equal to roughly 0-6, which also varies somewhat with 
the chemical state of the flow. The numerical differences are very small, 
however, when Le ~ 1.) 

The general result is seen to be very similar to the Couette flow result 


q hi bd {'H,—H,,.+ D(Le—1)(c,;—¢,,,.) + }Pru?', 

when allowance is made for the difference in definition of Nusselt number 
in the two cases and for the interpretation of 6. (The pu ratio dependence 
is not so clear cut, but this is hardly surprising since boundary-layer thickness 
is directly affected by the values of p and » whereas 6 in Couette flow is not.) 
The two solutions compare well in general shape when the flow is in chemical 
equilibrium and the wall is ‘cool’, for then ¢,,.~ 0 and H.>H 
H,+}3PrU?> H,,, in each case. An important difference arises when 
the flow is chemically frozen, however, as a result of ditferences in the 
interpretation of a fully catalytic wall. Fay & Riddell (in their later paper) 
remark that when the wall is catalytic, ‘‘the atom concentration will be 
reduced to its equilibrium value at the wall temperature’’. In particular, 
if 7’. is low enough, c,,,. would approach zero on this hypothesis, whereas 
it is clear from the present work that this can never occur since it would 
require an infinite value of the recombination coefficient I. When 
comparing q,., and q,,, values for Couette flow, the consequences 
of assuming c,,. = 0 in place of putting I’ = 1 are not too serious (in fact 
Guy Uwe fOr C,,. = 0 1s roughly 0-9 in place of 0-85 when [= 1). But 
Hirschfelder (1956) suggests that I’ for recombination of oxygen atoms 
on a metal surface is roughly equal to 2/7’, so that it is never greater than 
10-? in a practical case. For the examples quoted here, [ ~ 10-3 when 
T,,. = 1000° K and the lower wall is metallic, whence it can be seen from 
equation (57) that ¢,,. = c¢,,; = 1:0. Thus only a very small fraction of the 
dissociation energy is recovered and q,,,q,..~ 0:21, representing a 
considerable reduction in the frozen flow heat transfer rate. 


8. CONCLUSIONS 

It has been shown how the homogeneous chemical reaction rate 
exercises considerable control over the detailed properties of a layer of shear 
flow. When these reaction rates are fast compared with the rate of diffusion 
of the various chemical species in the mixture, local chemical equilibrium 
composition is closely approached and the nature of the surface reaction is 
is unimportant in the majority of the flow. In particular, energy transfer 
rate is effectively uninfluenced by the surface reaction rate. When the 
lower wall is sufficiently cool, as it may have to be in practice, the atom 
concentration falls to zero some distance above the surface and energy is 
finally transferred to the wall by thermal conduction alone. 

When homogeneous reaction rates become very slow compared with 
atom diffusion rates, the picture changes considerably. Energy transfer 
into the wall occurs through both possible mechanisms (even when the 


wall is cool), the thermal conduction contribution remaining sensibly 
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independent of surface reaction rate, whilst the atom diffusion contribution 
depends critically on the value of surface recombination coefficient. If 
the wall is a very effective catalyst for the recombination reaction, almost 
all the atoms which strike it will recombine and give up their dissociation 
energy to it directly. Thus a large ‘atom sink’ is set up at the wall, which 
requires large concentration gradients to provide the necessary rate of atom 
diffusion. For any given value of the ‘external’ concentration (c¢,,) an 
increase in the value of I’ will therefore result in lower atom concentrations 
at the wall, and hence in higher rates of final energy transfer. By the nature 
of the surface reaction, c,,. can never be zero there since, if it was zero, 
there would be no atoms present to take part in the reaction. 

It is suggested that even for quite high surface reaction rates (e.g. one 
atom in every 100-1000 atoms which strike the wall recombining) energy 
transfer occurs significantly more slowly in chemically frozen flow than 
in otherwise comparable circumstances with chemical equilibrium flow. 
It seems reasonable to suppose that the influence of the surface reaction 
spreads out into the shear layer as homogeneous reaction rates decrease. 


The present work is a revised version of a report written for the English 
Electric Co. Ltd., G.W. Division, Luton (Clarke 1957). The writer is 
particularly indebted to Dr M. McChesney, now with the B.T.H. Co., 
Rugby, for many valuable discussions on the problem. 
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SUMMARY 

This paper examines the flow characteristics of a body of small 
slope planing at high Froude number over a water surface. An 
equation is obtained relating the slope of the planing surface to 
an integral containing the pressure distribution on the planing 
surface. “The equation is expanded for large Froude number and 
a solution is obtained by an iteration process. At each stage of 
the iteration process the integral equation of ordinary thin aerofoil 
theory is solved. The pressure distribution on the planing surface 
is derived as a series in inverse powers of the Froude number F, 
as far as the F-* term. Computations are performed for the 
planing of a flat plate, a parabolic surface, and a suitable linear 
combination of these shapes which results in a flow without a 
splash at the leading edge. 


1. INTRODUCTION 

When a surface craft moves through water at low speeds it remains 
near its floating position relative to the undisturbed water surface, the 
lift being supplied mainly by hydrostatic forces. If the speed of the craft 
is increased through a certain critical value, the craft rises and moves over 
the water, with separation of the flow taking place at the trailing edge. 
The craft is then supported by hydrodynamic forces and 1s said to be planing. 

One of the features of a planing motion is the narrow jet of water, known 
as the splash, which ts thrust forward by the planing surface near its leading 
edge. Neglecting friction, the drag force on a planing craft is made up of 
two parts: (1) a splash drag due to the work done in creating the splash, 
and (2) a wave drag representing the energy carried downstream in the 
wave motion set up by the passage of the craft. If the craft, which is taken 
to be of length 2/, has a velocity U, waves fixed with respect to the craft 
are of wavelength A = 27F?/, where F = U(g/)~!? is the Froude number. 

This paper is concerned with the determination of the two-dimensional 
flow past planing surfaces which are inclined at small angles to the 
undisturbed water surface. A solution applicable at high Froude number 
is obtained for a planing surface of arbitrary shape. It is seen that the 
surface shape can be chosen to give a smooth flow at the leading edge, 
thereby eliminating the splash drag. Since the wave drag is small at high 
“roude number, the shape which eliminates the splash drag may provide 
good approximation to the optimum design for minimum drag. 

The solution in the infinite Froude number case has been given by 
Wagner (1932). The equations governing the flow are shown to be the 


— 


~ 
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same as in the flow past infinitely thin aerofoils. Except for the splash, 
which is assumed thin, the flow set up by the planing surface is identical 
to the flow in the lower half-plane of an unbounded fluid disturbed by an 
infinitely thin aerofoil of the same shape as the planing surface. The 
pressure distribution on the planing surface is the same as that on the lower 
side of the corresponding aerofoil and the lift is therefore half that of the 
aerofoil. An account of the method used by Wagner for obtaining the 
splash drag from the aerofoil solution is given in §3. Wagner examines 
the flow past a body of circular camber at various angles of incidence 
of the chord line. At zero incidence drag-free planing is achieved since 
the wave drag is zero at infinite Froude number and for a symmetrical 
surface there is no splash drag since the flow must be smoothly joined to 
the planing surface at the leading edge as well as the trailing edge. 

Lamb (1932, §§ 242-4) considers the two-dimensional flow due to the 
application of a pressure point to the surface of a stream. The planing of 
a body at arbitrary Froude number can be represented by a distribution 
of pressure points whose strength is proportional to the excess pressure 
on the body at each point. The boundary condition of zero normal velocity 
on the body then leads to an equation relating the slope of the planing surface 
to an integral involving the pressure distribution along the planing surface. 
Lamb considers two simple pressure distributions for which the integrals 
can be evaluated to give the shape of the planing surface. 

Two attacks on the integral equation involved in the inverse problem 
of finding the pressure distribution on the planing surface in terms of the 
given surface shape have been made by Maruo (1951) and Squire (1957). 
Maruo assumes a Fourier series expansion for the pressure distribution 


in the form 


1—cos 0 J . 
p = a4 ——— + } a, sinné, (1) 
sin @ wat 
where x = — /cos@ for 0 < 6 < z denotes the distance measured from the 


centre of the planing surface. An infinite series for the body slope results, 
comprising an infinite set of functions multiplying the coefficients a,,. 
From this set an orthogonal set is ccnstructed, which enables equations 
for the a,,’s to be derived. The example of the planing of a flat plate at 
general Froude number is computed, yielding the pressure distribution 
and the lift and drag coefficients as functions of the Froude number. 
A diagram is given showing the variation of the constants a, to a, with 
Froude number, but no other indication is given regarding how many 
constants have to be calculated for accurate results. Maruo’s work provides 
the solution to the planing problem of an arbitrary surface shape at general 
Froude number. However, the analysis necessary in obtaining and using 
the orthogonal functions described above is very complicated and the 
present paper is intended to provide a simple alternative to this method. 

Squire considers only the first four terms in a similar series for the 
pressure and calculates the coefficients to give a surface shape as nearly 
flat as possible by imposing mean curvature conditions. The results obtained 
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agree over a wide range of Froude numbers with those given by Maruo 
for the flat plate. 

A new approach is presented in §2 of this paper. This enables the 
solution for an arbitrary body shape to be easily computed, provided the 
Froude number is moderately high. It is seen, from comparison with the 
pressure distributions on a flat plate given by Maruo, that the approximate 
solution derived here provides accurate results for Froude numbers greater 
than 3. 

‘The equation relating the planing surface slope to an integral containing 
the pressure distribution on the planing surface is solved by an iteration 
process. ‘The integral is expanded as a series in inverse powers of the 
Froude number. An iteration process provides successive terms in the 
pressure distribution expansion in a similar series. At infinite Froude 
number, the solution to the integral equation is the thin aerofoil solution, 
vielding the first term in the pressure distribution expansion as an integral 
of the body slope. At each stage of the iteration procedure the same thin 
aerofoil theory integral equation has to be solved. The solution of this 
equation provides each new term in the expansion of the pressure 
distribution as an integral of the previous terms. ‘The expansion up to 
terms of order /* is obtained in §3 and the multiple integrals obtained 
in the solution are reduced to single integrals of the body slope, involving 
certain functions which are given in tables 1 and 2. The lift and drag 
coefficients are obtained in similar form. 

The large Froude number approximations for the planing of a flat plate 
and of a parabolic camber shape are given in §4. As already mentioned, 
a parabolic surface at zero incidence, planing at infinite Froude number, 
has smooth flow at the leading edge and therefore has no splash drag. 
An important result, that emerges from the study of the flat plate and 
parabolic camber flows at high Froude number, is that there always exists 
a linear combination of the two flows which represents a flow having no 
splash. In this high Froude number range, where the wave drag is small, 
elimination of the splash drag causes a considerable reduction in total drag. 
The planing characteristics for a body shape derived in this manner are 
discussed in § 4. 


2. LARGE FROUDE NUMBER APPROXIMATION 

A surface craft, having large span so that the motion can be assumed 
two-dimensional, is considered to be planing at a speed U over a water 
surface. ‘The water is considered to be incompressible and of infinite 
depth. ‘The planing surface is assumed to be inclined at a small angle to 
the horizontal at all points of its length. Cartesian coordinates are taken 
with the origin at the middle of the planing surface, the x-axis in the direction 
of motion, the y-axis vertically upwards and half the length of the wetted 
surface is taken as the unit of length. 

For an irrotational motion, a solution to Laplace’s equation is sought 
in the lower half-plane which satisfies the boundary conditions on the free 
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surface and on the body. The boundary conditions are linearized with 
respect to the perturbation velocities and are applied on y = 0, since the 
height of the water surface, 
y = n(x), (2) 
is of the same order as the perturbation velocities. 
The linearized boundary condition of zero normal flow across the water 
surface is 


dyn _ _ cd F 
_ ae (3) 


where UVd(x,y) is the perturbation velocity. On the planing surface 
v <1, dy/dx is prescribed as S(x), the slope of the surface. Bernoulli’s 
equation gives a second relation between 7 and ¢ to be applied on y = 0 
as 


( 


De a = + Ky + pd —— | 3 (4) 


CX 

where AK = F? and p is the density of water. The pressure p, measured 
in excess of atmospheric pressure, is zero on the free surface and is to be 
determined on the planing surface. The term p¢ in equation (4) represents 
Rayleigh’s artificial friction, which ensures no disturbance far ahead of the 
craft. The term can alternatively be regarded (according to Laplace 
transform theory) as the contribution for large time from ¢cd/ct, for the 
unsteady problem of a craft starting from rest. The required steady solution 
is obtained by eventually letting » tend to zero through positive values. 

A potential which satisfies the free surface conditions is given by Lamb, 
the planing surface being represented by a distribution of pressure points, 
and may be written as 


, 1 . ur ty yy 
bs Raa » -Pleye | sana * dé, (; 
~~ we Jo 4-K+tp 


where p(€) = pU* P(E) is the excess pressure on the body. ‘The boundary 


41 
— 


condition (3) on the planing surface leads to the integral equation for the 
pressure distribution 

: 1/7 ., dé 1K / : "x gia(r— §) ; 
S(x)=-} P(€)—— + lim #A7— | P(E) | ——— dx dé, (6) 
1 ee >t) 7 I /0 4—K-4 1h 


where the Cauchy principal value is to be taken in the first integration on 





the right-hand side. 

The inner integral in the last term of equation (6) will now be transformed 
into the exponential integral by putting 

v—-€ =A, ta(x—K+ip) = X. (7) 

‘The integral with respect to Y parallel to the imaginary axis is deformed 
into an integral parallel to the real axis, as shown in figure 1, there being 
no contribution from paths at infinity in the negative real half-plane. 

On letting » —-0 through positive values, 


in(r £) 


Dg ee paiaB gX ; 
lim \ Ha. dx = —e (20iH(—a)+ | <x); (8) 


4-0) - tL 
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where H(x) is the Heaviside unit function. The residue term —2zie'“* 
is included when a < 0 since the contour is then deformed over the pole 
at XY = (0. The expansion for small K of the integral on the right-hand side 
of equation (8) is given by Jahnke & Emde (1945), enabling equation (6) 
to be written 
: (" P(E) ee 
. 1 Saal” | 

K /} : 

— | P(€)[~7—4dasgn(x—€)+(y—1)AK(x-€)+ 

7 1 

+K(x—€)log(K x-€)+...] d€&, (9) 


where y = 0-5772 is Euler’s constant. This equation is solved by an 


iteration process on the small parameter A, yielding the pressure P(€) as 


S(x) = 


a series in K. 


a 
Spake ancienlaglh chase v0 


(a) (b) | 


Figure 1. Integral paths for the exponential integral in equation (6), travelling from 
the point P = —1a(K-—i) parallel to the imaginary axis, are deformed into 
paths parallel to the real axis as indicated in (a) fora > O and in (4) fora < 0. 
The contour is deformed over the pole at X = 0 in (8). 


In the infinite Froude number case the integral equation (9) for the 
pressure distribution on the body reduces to 


“1 
Sy=-{ Pe) = 


1 &—x 





(10) 


This equation is that of thin aerofoil theory and is discussed, for example, 
by ‘Tricomi (1957). The solution to equation (10) gives the first term, 
P,(€), in the expansion of the pressure, as 





‘ 1] 1+é€\ (2 S(t) 1—t' 
A) cod Sa: poked P'S 
o(€) I) ~G)« (11) 


Equation (10) admits another solution of similar form which is rejected 
on account of the Kutta condition that the singularity has to be at the leading 
edge. Using the solution for P,(€), equation (9) when approximated to 
the first order in K becomes 


S(x)-K | Py(é)[1—-4sen(x-Qy ae = 2 /° Pey&, 12) 
J -1 7 


Ts a3 €-x 
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where the known functions have been collected on the left-hand side. 
This integral equation has the same form as the aerofoil integral equation. 
\ solution for the second approximation P)(é)+ AP,(&) follows, where 
P,(€) is given by 


— 1 1+é\/3 a — 1—t\ dt ‘ 
pie) == /(4)| [| Pa(ry(—ssen(e—myar |, [(F=1) ©. (13) 


Equation (9) may now be approximated to the second order in A using the 
above solution for P)(é)+ A P,(€) and the third approximation 


P(E) + K Py(E) + K* P(E) 


to the pressure distribution follows in like manner. P,(€) is given by 
1 1+€\ /} be 1—t\ dt 
P,(€) = - - O(r, t) dr — : 14 
ale) 7 at L ee yer! Maan oa 


O(r,t) = Py(r)(1— } sgn(t—r)) + 


Po(r)[(y — 1)(t—1) + (t—r)log(K t—r )). 








where 


The expressions (11), (13) and (14) provide the expansion of the pressure 
distribution on the body up to terms of the order of F~“ in the form of 
integrals involving previous terms in the series. In the following section 
these integrals are transformed by successive substitutions to ones con- 
taining the given function S(t) only. The resulting multiple integrals are 
reduced to single integrals. 


3. SIMPLIFICATION OF THE SOLUTION: LIFT AND DRAG 

In the expressions for the terms in the expansion of the pressure distribu- 

tion on the body given in the previous section, the distance integrations 

along the body can be replaced by angular integrations by substitutions 
of the form 

€= —cosd for0<d<z. (1 


wn 
— 


Equation (11) for the first term in the pressure distribution expansion may 
then be rewritten 


g 6 
P,( = cos ¢) = A, tan ld oh. d 


<j [” te esl —__—_—_. ¢ 
ind i: S(—co a sie 
( S(—cos@) dé. 


- 0 


Ay = 


where 1 


Substituting this expression for P)(—cos¢) in equation (13) for the second 
term in the expansion leads to 


“a 7 1 4 
P,(—cos¢) = A, tan}¢+B, T(¢)—- bs | a 


XK 
7 | LJ 9 cos8—cosx% 


™ §(—cosy)sin?y d) is 
«(f (—cosy)sin per) ae | (17) 


cos §— cosy 





! 








472 E. Cumberbatch 


\ ‘re ? lied 
where 1, E \7B, ae | S( —cos #)sin fa] T(@) dé, 
TO 
9 “7 
B, = = | S(—cos4#)(1+cos4@) dé, 


| log tan 36 dé. 


0 


T(¢) = - 


Values of 7(¢) are given in table 1. This form of P,(—cos¢) is convenient 
for some very simple S(—cos¢) since the formula 


cos né zsinnd ' 
| oe dé = -_—_ (18) 
|) cos9—cosd@ sind 


enables the multiple integral in (17) to be evaluated directly. For general 





d, T(h) Zh) b T(d) | Z(dh) | 
0 0 0 | 26 0°-3579 0-1369 
1 0-0319 0-0355 | 30 | 03887 | 01345 | 
») 0-036] 0-0541 3 04165 | 01305 | 
3 0)-0774 0-0678 38 0-4415 0-1251 
4 0-968 0-0788 4? 0-4640 0-1185 
5 )-1148 0-O878 46 0-4841 0-1111 
6 00-1316 | 0-0954 50 0-5021 | 0-1030 
7 0-1476 0-1019 54 05181 | 00942 
R ()-1627 Q-1075 38 00-5321 0-O849 
ov) 0:1772 0-1124 62 0-5443 0-0751 

10 0-1910 0-1166 66 05548 0-0650 

12 Q-2171 ()-1234 70 0-5635 0)-0546 

14 0-2412 (): 1284 74 05706 0-0440 

16 0):2637 0-1321 78 0-5761 0-033 

18 0): 2848 ()-1347 aS 0-5800 0-0222 

0) 0-3047 0-1364 86 0-5823 0-O111 

2? ():3234 0-1372 90 05831 0 


‘able 1. Values of the functions T(4) and Z(d) defined in equations (17) and (21), 
with the argument d measured in degrees. For the range 90 < 4 < 180°, 


the formulae 7(d) T(180—¢4), Z(d) Z(180—¢) may be used. The 
second of these relations follows from the result | log? tan}é dé = }7°. 


body slope the multiple integral in (17) is reduced to a single y-integration 
suitable for numerical computation by first performing the @- and w- 
integrations. ‘Che changing of the order of integration requires caution 


because of the double pole in the integrand when 6 = y = J. The identity 


sin 6 * S(—cosy)sin? y 
| —_—_—— } a), SLRs Ay ates : dy dé 
9 cos 0 — cosy 90 cosé6—cosy 


sin @ 7 (S(—cosy)— S(—cosd))sin? y 
teeta | ( y) ) : dy dé + 
0 


] cos §—cosy 











'g coS4— cosas 


ad Sl fd w Sl 2.4, ‘ 
+ S(—cosy) | = I snl ay | dd (19) 





9 cos @—cos is 9 cos 6—cosy 
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overcomes this difficulty since the order of integration in the first term 
on the right-hand side may now be changed as it does not contain the 
double pole. The second term in (19) may be evaluated directly since the 
integrations do not involve the body slope. Performing the @- and ¢- 
integrations and making the substitutions f¢ = tan}y(tan}y)-! and 

u = tan }d(tan }y)~', equation (17) is reduced to the form 

P,(—cos¢) = A,tan}¢+ B, T(¢)+ 
} = tan® 3¢ | u Fo See a¢ = ma) (u* + tan? 1d)-? — 
i Jo tan? 36+? 

a s(3 tan? }f—1 
u* tan? 3¢+1 


oe tan? }4¢+1) ‘| du, (20) 




















where F(x) = | oat dt. 
ler= ] 
Values of F(x) are given in table 2. 
| | 
x F(x) — F(x) | x | F(x) — F(x) 
| 

0 0 | 9 | 016 | 0-4562 1-4504 
0-001 | 00079 | 00635 | O18 | 04927 15151 | 
0-002 | 00144 | 01061 | 020 | 05272 1-5724 
0-004 | 00261 |  0-1741 | 0-22 | = 05598 1-6234 
0-006 | 0:0367 | 0:2304 — | (0:24 | 5909 1:6691 
0-008 0:0466 | 02798 | 0:26 |  0-6205 1:7101 
0-010 | 00561 | 03242 | 0-28 | 06487 | 11-7527 
0-015 | 00780 | 0-4206 | 0-30 | 0:6758 |  1:7806 
0-02 0-0983 0-5026 0-35 | 0-7386 | 1-8514 
0-025 | 04172 | 05747 | 040 | 60-7957 |  1-9075 
0-030 | 01352 | 06394 | 045 | 08480 | — 1-9523 
0-035 | ()-1524 ()-6982 0-50 0-8961 | 1-9881 
0-040 | 01688 | 0°7522 | 0-55 | 0-9406 | 2-0168 
0-045 | 0-1847 | 08022 | 0-60 | 0-9819 | 2-0398 
0-050 0-1999 | 0-8488 0-65 1:0205 2-0597 
0:06 | 02290 09-9332 | 0-70 | 1:0566 | 2-0721 
0:07 | 0-2565 10084 =| 0-75 | 10905 | 2-0831 
0-08 | 02825 1:0759 | 0-80 | 11224 =| = 2-0912 
009 =| 0:3074 | 21371 | 085 | 11526 |) 2-0970 
0-10 0°3311 | 1-1931 0-90 1:1811 | 2-1008 | 
0-12 | 03759 | 1:2919 | 0:95 | 1-2081 | 21023 | 
O14 | 04174 | 13767 | 100 | 1-2337 | 2-1036 





Table 2. F,(x), F(x) are defined in equations (20) and (21). The relations 


F\(1) = 147, F.1) = —}3a | T(d) dd, provide a check on the accuracy of 


tables 1 and 2. 

P,(—cos¢) may be similarly reduced to single integrals involving the 
body slope. The logarithmic term in equation (14) may be dealt with by 
integrating equation (10) from —1 to y with respect to x and then from 
~1lto r with respect to y. Using techniques similar to those applied in 
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reducing P,(—cos¢), P;(—cos¢) may be expressed as 
P,(—cos¢) = A,tan}¢+B, T(d)+C,sind — (24,/7)sin d log tan $4 — 
— (B,/7*)sin d log? tan 36 + 

*. S(—cos6)sin® 6 aw) 24 )— 


“0 


(198, cosd— | 


i6 . ‘ ‘1 , 
— — singtan® 3d | (u—u?)F,(u) >» 
7 /0 


(tan? id—u?\, , 
x | s( gra) + tan? $46) — 





tan? t64+u? 
oe ee! ae 
-s(* — j)(wPtan® 1d +1) | du (21) 


u* tan? 36+ 1 


where 
9 - has 
1, = — 241 4 72 + [" 5(—-cos0)¥ (0) ad, 
- Zz 0 


2B 2 We 
B, = -—2A,- — - = | S(—cos #)sin? 6 log tan $6 dé, 
7 7 J0 


C, = —47LB,+ | S(—cos6)sin6Z(6) dé, 


/0 


(2y- 1+ 2log }K+ 4 |, 700) a), 


oiies 7 sa 
ati \ “7.0 
Y (6) = (2/7*)sin® 6 log* tan 36+ L sin? 6— } sin 24 Z(6), 
, SE 6 
Z(6) = — | log*tan}iy dy— —, 
am 1 -/ ‘ : = 
Jo 27 
; log? t 
Fi(x)= | =~ dt. 
2_ ] 
; 


Values of the functions Z(#) and F,(x) are given in tables 1 and 2. 


The lift coefficient, 


Cy = P(—cos ¢)sin ¢ dd, 
0 
may now be obtained to the second order in A as 
: B C, B. B 
. ee Se Se Bek ll AS a SS ee 
C= <r | =( 42 t+ St)s A 


log? tan 10 )sin* cos 4 a | K?. (22) 


by 1 
_ oe ‘os = het 
] 0 ( ‘ (Gj 7 


Similarly the drag coefficient, 


C, = | P(—cosd)S(—cosd)sind dd, 
to the second order in A is 
t[ A? + (2A, A, + 40 B?)K + (24, A. + A} + $7B, B,)K?). (23) 


C cae 
The splash drag is determined by considering the singular term A tan }¢ 
Wagner interprets this term, 


in the pressure distribution on the body. 
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which corresponds to a term with a z!? behaviour in the complex potential, 
as describing the reversal of the fluid in the narrow jet of width 6 forming 
the splash and he obtains the relation 6 = }74?. The work done to produce 
the splash may then be calculated in terms of 4 to give the splash drag 
coefhicient as 


C, = #A*. (2+) 
To the second order in A, therefore, 
Cy, = 7[A2+2A, A, K+ (24) A. + A7)K?]. (25) 
The shape of the free surface is 
n(x) = | . P(€)(1 —sgn(x—€))sin K(x — &) d& + 
J —I 
i; — 
+ - P(€)[{ 37 sgn(x — €&) — S,(K(x— &))}sin A(x — &) — 


3 

~C,(K(x—€))cos K(x—£)] dé, (26) 
where S; and C, are the sine and cosine integrals. The wave drag is derived 
from the form of the waves far downstream. The first term on the right-hand 
side of the equation for »(xv) represents the wave train downstream, while 
the second term represents a local disturbance. The wave drag coefficient 
is deduced to be 


“1 2 “J 2 
Cy = | ( | P(€)cos KE dé +( | P(£€)sin KE ae) | 
= 1 . 1 
= }7?B? K+ 40°B, B, K?, (27) 
expanded to the second order in A. It is noticed that the large Froude 
number approximations for the splash and wave drag given by (25) and (27) 
add up to give the total drag (23). The relation 
Cy = KCl (28) 


valid to the second order in A, may be obtained from equations (22) and (27). 


4. LARGE FROUDE NUMBER PLANING FOR FLAT 
AND PARABOLIC SHAPES 

In this section the large Froude number expansion for the pressure 
distribution up to terms of order K? is derived for planing surfaces of flat 
and parabolic shape. For the planing of a flat plate the body slope is given 
by S(—cos@) = a, where x is the inclination of the plate. Substitution of 
this value of S(—cos @) in the expressions (16), (20) and (21) for the pressure 
distribution in the previous section gives 


(1/x)P(—cos¢) = A, tan }¢6+B, T(d)+C,sin¢d— 


— Fisk. ; 
- (= K -2(1 ~ =) K+) sin blog tan Io+ 


+ (K?/z*)sin 6(2 + cos d)log* tan 44 + 
+1K*sin2¢—47K7(1+2cosd)Z(d), (29) 


2H2 
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where 

: eee : ; 1 Les \ : 
A, l—(7+2/7)K+ }K*log }K+(7?+6+ = +hy+ - | ad! K* 
Bb, ~-2K+2(7+4/7)K?, 


Yibo 

~ 
— 

aa 

= 
— 

~ 

~ 

= 

x 

“———" 

to 


bol 


Ph, 
/ 


x 














° Po { 


Figure 2. Distributions of the pressure coefficient P/« on a flat plate at inclination x, 


planing at Froude numbers corresponding to K = 0,0:05 and 0-1. The 
lift coefficients in these three cases are C= 7a, 0-827« and 0-737, respectively. 
The leading and trailing edges are at x = 1 and x = —1. 


The pressure distributions on a flat plate planing at Froude numbers 
corresponding to K = 0, 0-05 and 0-1, are drawn in figure 2. Comparison 
with the exact pressure profiles on a flat plate obtained by Maruo show that 
the large Froude number expansion up to terms of order K®? is adequate 
for K less than about 0-1. At K = 0-1, which corresponds to a Froude 
number of about 3, there is a 6°, discrepancy. 

For a planing surface of parabolic camber the slope of the surface is 
taken to be S(—cos@) = 8 cos#. The pressure distribution expansion to 
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terms of order AK? for this shape is 
(1/8)P(—cos¢) = A, tan3éd+B, T(¢)+C, sind—- 


K K*\ . ss 
- (= — —~ }sin2¢logtan }¢+ K*sin¢ log tan 34 4 
oe 





r9 9 


oe fate a ) Se 
1 sin? ¢)sin ¢ log? tan 36 — PD sin? d — 


— = K*cos¢Z(h), (30) 





2 
where 
&, = —4tnK+(4n7+ 5) 2. 
B, = -K+(74 =) 
37 
C,=1- K +iK?lopiK+ ives ! . - ! (" T(@) dé) K?. 
A 7 s Pea 2 + 37? TJo 


At infinite Froude number the singular term A, tan }¢ in the pressure 
distribution disappears. This represents a flow which is smooth at the 
leading edge and the motion is achieved with zero splash drag. ‘The pressure 
distribution in this case is elliptical. 

From the pressure distributions for the flat plate and parabola, it can 
be seen how a flow without a splash drag may be obtained at any preassigned 
high Froude number. If the planing surface has a slope which is a linear 
combination of the slopes S(—cos@) =a and S(—cos@) = Bcos@, so 
chosen to eliminate the singular term tan $¢ in the resulting pressure 
distribution, then a splash-free flow will result. For the Froude number 
range F > 3, the pressure distribution on a planing surface chosen in this 
way approximates to the elliptical distribution on the parabolic shape 
planing at infinite Froude number. The required combination is obtained 
by choosing the constants « and f to satisfy the relation 


where A, and 4,, are defined in equations (29) and (30). It is to be 
emphasized that the splash-free flow obtained by satisfying this relation 
is attained at only one particular Froude number. As f determines the 
shape of the planing surface, the above relation is satisfied for any given 
parabolic shape by trimming the craft to give the required angle of 
incidence z. 

The lift coefficient for a combination of the flat plate and parabola is 
derived by substituting S(—cos#) = «+f cos @ in (22), and is given by 

C, = an(1—4-41K + K* log }K+19-9K?) + 

+43Br(1—3-99K + 3K? log3K+17-4K?). (32) 


L 


This represents the sum of the lift coefficients for the flat and parabolic 
shapes (the pressure distribution for the combination is derived by adding 
the pressure distributions in (29) and (30)). Equations (31) and (32) may 


/ 
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alternatively be used to express x and f in terms of C, and K (correct to 
terms of order K?) in the form 

z= KC,, (33) 

B = (2/7)C,(1+0-85.K — 3K?log }K —0-2K?). (34) 
These equations give the necessary incidence x and shape (specified by 8) 
to obtain a given lift coefficient for planing without splash at any particular 
Froude number F > 3. For this range of Froude numbers, corresponding 
to K < 0-1, the A? terms in (34) are less than 2°, of 8, and both x and 8 
may be taken as linear in K. 

The shape of such a splash-free planing surface for a Froude number 
corresponding to K = 0-05 is shown in figure 3, together with the local 
wave profile. Since the lift coefficient determines the magnification in the 
vertical direction it is convenient to plot y(C,)~! instead of y. The shape 
of the free surface is obtained from equation (26) by numerical integration 
using the large Froude number approximation for the pressure distribution. 
The drag coefficient for the splash-free craft planing at K=0-05 is derived 





re} 30 j os x 
-2 
Figure 3. Splash-free craft planing at K = 0-05 shown with local free surface shape. 
The craft lies between x = —1 and x = +1. Vertical scale is enlarged as 
indicated. 


from (28) as C, = 0-05C}. By way of comparison, the drag coefficient 
for a flat plate planing at the same Froude number is C, = 0-39C?. This 
drag is obtained by adding the wave drag (28) to the splash drag C, = 7A3, 
(24), where A, is defined in terms of x by (29) and may be expressed in 
terms of C, by using (32) with 8 = 0. A considerable reduction in drag 
is thus effected by eliminating the splash, so that the shape of such a 
splash-free craft may provide a good approximation to the optimum design 
for minimum drag. 


The writer wishes to express his thanks to Professor M. J. Lighthill 
and to Dr R. F. Chisnell for many helpful discussions on this work, which 
was carried out while the writer held a Department of Scientific and 
Industrial Research maintenance grant. 
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The mixing of ground water and sea water in 
permeable subsoils 
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SUMMARY 

The subterranean mixing in permeable media of sea water 
and ground water is studied. ‘The model for this mixing process 
which was suggested by C. K. Wentworth is adopted, but is soon 
discarded in favour of a more tractable formulation whose 
equivalence to the original model is established. The analysis 
is carried to the point where the determination of the salinity 
distribution of the ground water in a given subsoil requires only 
the solution of an elementary linear ordinary differential equation. 


1. INTRODUCTION 
The distribution of water in permeable islands, e.g. the Hawaiian group, 
has the general configuration shown in figure 1. The depth of the lens-shaped 
region of fresh water (usually called the Guyben—Hertzberg lens) is nearly 








y Y 


Figure 1. The distribution of ground water in a permeable island. The permeable 
material lies below AOB; the fresh water lies in the region between EO and 
DO; the salt water lies below DOB. A description of the salinity distribution 
near DO is the objective of this investigation. The arrows indicate qualita- 
tively the velocity distribution associated with the ‘ steady’ rainfall and 
run-off to the sea as described in § 4. 
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proportional to the one-half power of the distance from shore; its size is 
determined by the net water supply. However, the configuration is not 
stationary since there is an oscillatory vertical motion of the fluid in response 
to the tidal pressure excitation along OB. This motion implies, in the 
neighbourhood of the interface between the fresh water and the salt water, 
an invasion by the salt water of the cellular medium assigned to the fresh 
water, and vice versa; this leads to a dispersion of the material, so that the 
interface becomes diffuse. The purpose of this paper is to investigate 
quantitatively the structure of the transition zone separating the salt water 
from the fresh water. 

In §2 a physical model for the mixing process (due to C. K. Wentworth) 
is introduced; a continuum model which is thought to be equivalent is 
also suggested. The equivalence is established in § 3, by comparing solutions 
to certain preliminary problems. In §4a more realistic continuum formula- 
tion of the problem is given, and a solution is obtained; from this a more 
convenient equivalent mathematical model is established. In §5 this last 
model is used to solve what I believe to be a satisfactory formulation of the 


problem of the structure of the mixing zone. 


2. "THE MIXING PROCESS 
We shall adopt a mechanism due to C. K. Wentworth to explain the 
salinity distributions associated with the problems of interest. Imagine 
the porous structure to be a homogeneous array of communicating cells; 
a one-dimensional array like that of figure 2 suffices for our present purpose. 


Figure 2. One-dimensional array of cells 


+ 


Denote by S,,, the salinity of the fluid in cell n at time ¢,,, and define a 
velocity W(t) = W pdf, where M is the upward mass flow of fluid per 
unit time in the array of figure 2, 4 the cross-sectional area of the array, 
p the fluid density, and f the porosity (i.e. the fraction of the volume that 
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can be occupied by fluid). According to this definition, W is the speed 
of a free surface advancing through the porous medium ahead of such a 
mass flow M. The cell height is denoted by A. S,,,, may be taken to 
be either the volume fraction or the mass fraction of salt in the solution. 
During the time interval f¢,, to ¢,,.,, a volume of fluid (t,,.,—t,,) WA 
is convected from cells n—1 into cells m. We assume this process takes 
place without diffusion and that all the fluid in a cell m mixes thoroughly 
at time ¢,,.,. We choose the time interval sufficiently small so that no 
fluid from a cell n—1 passes into a cell n+1, Le. t,,.,—t,, < A/W. 
Assuming that W’ is positive in the time t¢,,.,—t,,, cell m loses salinity 


S',.,(1—a) and gains salinity S,,,_,(1—a), so that 
Snctn— Sun = (1—@)[S,,.n-1—Sn.n}> (2.1) 
where l—a= W (t,,.,—-t,,)/A. 


When JV is negative during the time interval, S,, ,_, must be replaced 

by S,,,.,- To include both cases, the difference equation governing the 
Y Oman 1 
S',., Must be expressed as 
: @ (1 —a) , are ree ss sat 
6144 ~Sen 7 [((Ww+ W)S,,,1-2WS,,,+(W -—W)S,,, <1). 
(2 

Another model which has the same plausibility and which gives the same 
prediction, as we shall see in $3, can be constructed by the following 





2) 


conventional limiting process. Divide each side of (2.2) by t,,—t,,, 
and let this ditference tend to zero. We obtain 
S’(t)= {(W +W)S,_,+(W —W)S, .,-2 W'S,}/2A. (2.3) 


Here S(t) denotes the salinity in cell ” at time t, and the prime denotes 
differentiation with regard to t. The corresponding formal limiting process 
wherein \-—- 0 would give a model in which the dispersive process had 
been eliminated. Consequently, we shall merely postulate a continuum 
model which ‘resembles’ (2.3) and establish its ‘equivalence’ with (2.3) 
by demonstrating that its interesting consequences are the same as those 


of (2.3). This continuum model, which is found by associating 
S,.,—-2S,+5S,_, with A’S,,, S,.,-—S,_, with 2AS,, and nA with y, 1s 
S,= |W\AS,,— WS,. . 2 


Here S(t, y) is the salinity at time ¢ and coordinate y, and the subscripts ¢ 
and y denote partial differentiation with regard to those variables. 
Equations (2.2), (2.3) and (2.4) are each solved in $3 for a given group of 
problems, and the equivalence of ‘cir predictions is established. 

We here briefly discuss restrictions on the size of a. For a subsoil of 
complicated geometry this parameter must be determined experimentally, 
but for some special geometries we can deduce certain information about 
its value. Ifthe cell array is composed, as in figure 3, of a tube with arbitrary 
vertical subdivisions, and if the time interval we adopt is the maximum 
possible for (2.2) to be valid, then a = }, since the flow at the very low 
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Reynolds numbers experienced in these phenomena would be paraboloidal 
in profile. Also, the volume of a paraboloid between the vertex and a section 
normal to the axis is one-half of that of the enclosing cylinder. If we use 
a shorter time interval, a is larger than }. If the cells are similar to those 
of figure 2, the flow from cell n—1 is likely to penetrate cell m in a 
comparatively slender filament, so again a will be larger than $. On the 
basis of such observations, it appears reasonable to expect that a> }. 
This is not a crucial point, however, since the determination of the effective 
cell size will be at least as important in determining the ‘diffusion rate’ 
of the phenomena. 


| 
| 
| 
' 


Figure 3. Tube with conceptual subdivision. 


3. ‘THE SOLUTIONS OF THE DIFFUSION EQUATIONS 
The treatment of equation (2.2) with — x < n< «, and with S,, 
given, 1s most readily carried out when we introduce a generating function 
gen(z)= > S,,.2", 
where the expansion is a Laurent expansion. To solve (2.2) we multiply 
it by s” and sum over n. The result is 


W+wWw W-wU 
Em+i(2) = | 1-(d-ay(i-2) So" =({i-a}i=—s ) Sar emt) 


If m, is the number of time intervals during which W > 0, and m, = m—m,, 
the solution of (3.1) is 


2mn(2) = Zo(2) [1-1 —a)(1—s)) [1 - (1 -—a)(l-—s7))™. (3.2) 


Since S,,,, 18 the coefficient of the Laurent expansion of g,,(), 
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The contour encloses the origin and passes to the left of the point z = 1. 
For all the problems of interest here, m > 1; thus, an asymptotic evaluation 
of the integral (3.3) provides all the useful information. The solution 
depends on the initial conditions, 1.e. on the values of Sy, No generality 
is lost if we take the initial distribution of salinity to be the step function 
which is unity for 2 > 0 and zero for n < 0, 1.e. g9(z) = (1—z)7 

With g.(z) = (1-—2)77, 


se |" (1—2)-lexp{ —in6 + m, log[1 — (1-—a)(1—z)]+ 


S 


+m,log[1—(1—a)(1—2-")]} dé, (3.4) 


where 2(@) = e’’, and the path is indented to pass above the origin in the 
#-plane. The saddle point of the exponent in (3.4) lies at 

4 = 6, ~ —i[n—(1—a)(m,—m,)]/a(1—a)(m, +m.) 
when 

n—(1—a)(m,—my,) < m. 

Since 4) < 1 for such n, the second derivative of the exponent is closely 
approximated by —a(l—a)m. Thus, the integral which asymptotically 
approximates (3.4) is, with A = a(1—a)m, 


S 
j ~ 


1 
mn ym 
1 


| [8 +(8—4)]}° 1exp{4A[#? — (0 — 9)? ]} d(A — 4%) 


erfc| — [n—(1—a)(m, — m,)]/[2a(1—a)(m, + m,)]*? 


i - -/ 12 
7 berfed —| y— [ | W(r) dr |/| 29 | maay' L (3.5) 
l J0 / /0 J 


The final equality is obtained when we use the coordinate definitions which 
were introduced with (2.4). This result states that the mid-point of the 
transition zone translates with velocity I(t), and that its breadth is 


aL 12 
proportional to I W(r) ar | : 
/0 
The generating function G(t,z) for the S,,(t) in (2.3) has the form 


G(i,z) = y S,,(t)z" = G(0,z)exp{ —2(1-—2z)-B(1-—2z7)}, (3.6) 


where x = (2A)-1 |" (W|+W) dt and 8 = (24) | (W|-W) dt. From 
/0 0 


this, the integral representation for S, (t) is found to be precisely the limiting 
form of the right side of (3.4) when a —> 1, i.e 


: 1 ; “ 
S,(t) = = | ‘ (1 —z)-lexp{—«(1—2z)—8(1—27-!)—in6} dé. (3.7) 
A discussion of the saddle point location and the second derivative of 
the exponent for m < x + shows that (3.7) leads again to (3.5). For large n, 
this could not be because the signal speed of the model of (2.2) is finite and 
that of (2.6) is not. 
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We now denote the solution of (2.4) by S(t, y) and its Fourier transform 
by a(t,7), 1.€. 


a(t,y) = | S(t, v)exp[—my] dv. (3.8) 


By the use of conventional techniques on (2.4) we obtain 
a(t,n) = [—iWy-— W Ay? ]o(t, ). (3.9) 
The solution, with S(0,y) equal to the unit step function, is 
a = o(0, 7 )exp{ — (% + B)A*y? —7(x— B)An}, 


where « and f are defined following (3.6), and o(0,7) = (m™)7. Thus, 
, ne —— ; a 
S(t,y) = 5- | (in) exp[ — (% + B)A2n? —7(% — B)yAX+ my] dy. (3.10) 


This time equation (3.5) rigorously defines the function implied by (3.10). 
The foregoing comparisons could be extended to more general initial 
considerations, but no further justification for the use of the continuum 


model seems necessary. 


4, MIXING WITH A SPACIALLY VARYING VELOCITY FIELD 
‘The Guyben—Hertzberg lens of fresh water is maintained in the presence 
of rainfall and the accompanying drainage of fluid to the sea. Thus, the 
velocity distribution above the salt—fresh water interface must resemble 
that of figure 1. Because of the diffuse character of the interface, some of 
the run-off must be salty, whereas the vertical intake is fresh water. If there 
were no motion other than the tidal fluctuations below the interface, the 
concentration at any point (for any initial solute distribution in dynamic 
equilibrium) would diminish with time. However, such a diminished 
concentration could not continue to be a hydrodynamic equilibrium 
configuration, and more of the denser salt water would have to be supplied 
by the sea to restore equilibrium. ‘Thus, the only acceptable velocity 
distribution both above and below the interface is depicted in figure 1. 
It is easier to generalize the ditfusion equation (2.4) to deal with an appro- 
priate one-space variable problem and then tackle the problem with the 
flow field of figure 1, than to deal immediately with the latter problem. 
The appropriate one-dimensional problem is one in which the vertical 
velocity is given by 
We wv coswt—eyv. (4.1) 
‘The corresponding physical situation is one in which this velocity distribution 
occurs 1n a vertical column of cells and a flux of fluid emerges laterally from 
these cells with a mass conserving distribution. With this velocity 
distribution, (2.4) becomes 


>; (eV —vcoswt)S§ y+ EV —UCOS wt AS, ‘ (4.2 ) 


Because of the variable coefficients and the non-analvtic character of the 


coeticient of S,,, rigorous solutions to this differential equation would 
be very difficult to obtain. We must simplify it, and the motivation for the 
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simplification can be better appreciated if we first anticipate some physical 
features of the phenomenon. With the ey term absent the solution of (4.1) 
is given by (3.5) and is 

S(t, vy) = Serfe{ — Ly — (v/w)sin wt]/[4t7— + 2P(t)]}* (4.3) 


where P(t) has period w and average zero. Hence .S has an error function 
distribution with a breadth proportional to t'? which translates at speed 
v cos wt and about which fluctuations occur at frequency w. We can expect 
that the principal effect of the ey velocity contribution will be to squeeze 
the mixing zone into a narrower configuration, and that the other features 
will remain qualitatively unaltered. 

In order to simplify (4.2), we note that the term jey—wcoswt| is, for 
all ey < v, very well approximated by wvcoswt|, except when t is very close 
to }(2n+1)z. Since almost no mixing occurs during such time intervals, 
the final term of (4.2) can be replaced by A'vcoswt|S,,,, with negligible 
loss of accuracy. When we introduce the dimensionless quantities + = wt, 
v = y(e/vA)!?, x? = eA/v, B® = w*A/ev, equation (4.1) becomes 


BS, = (ax—cos7)S,+a cos7|S,,. (4.4) 


We also write € = x—8-1(1+N~*)-'sin[7+arctan N-!] with N = B/a. 
The transformation to € allows for the oscillating translation of the salinity 
distribution. Equation (4.4) becomes 


NS, = €S,+ \cos7/S,,. (4.5) 


Typical estimates of the parameters are 8 = 107, x = 10~*, N = 10* (hence, 
terms of order \~! in the definition of € can be dropped with no important 
loss of accuracy). Since we are primarily interested in the salinity distribu- 
tion after a long time with the boundary conditions S > 1 as > x, S>0 
as €-» — %, we again use the initial condition at t= 0 that S = 0 for 
negative € and S = 1 for positive €. 

If we define the Fourier transform 


S(o,7) = | e-™* S(E,7) dé (4.6) 


with o = mm, (4.5) becomes 


NS_ = —(oS), + cos7 oS. (4.7) 
Furthermore, with z = logo and 
u=7—2/N, o=7+2/N, (4.8) 
(+.7) can be written 
2NS, + [1—|cos }(u+v)|e?-]S = 0. (4.9) 


The initial condition at ts = 0 becomes S = 1/o when u+v = 0. Equation 
(4.9) is a first-order linear equation which is solved by using an integrating 
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factor. That solution which is consistent with the foregoing initial 
conditions and in which the substitutions (4.8) have been used is 





S(o, 7) = 0 exp[o?/H(7)], (4.10) 
where 
H(r) = 7—1(1 —e-27' 4) + N-1[P(7) —e-?"5 P(0)]. (+.11) 
The function P(r) is a periodic function whose Fourier series is 
2 cos(2n7— ¢,,) 
Nol = FS (—1)*-1 " : 4.12 
re wad <4) 7(4n? — ner + N-2)12 ( ; 
where ¢, = arctan(n). Note that P(r) < > 2/(37n*) < 1. 
Equation (4.10) can be inverted to give a 7). The result is 
S(é,7) = 4(1+erf{4é[A(r)]'*}) (4.13) 


Recalling that P(r) < 1 and that a typical value of .V is 10', it is clear 
that after a large time 
S(é,7) ~ 3[1 +erf{(47)!78}] 
}[1 +erf{(47)!2[x— 8! sin r]}. (4.14) 
It can now be seen that when NV is large and when only the results for 
large + are wanted, (4.5) can be replaced by one in which the ‘ diffusion 
coefficient’ cosz is replaced by its average value 27. Note, however, 
that this asymptotic behaviour is reached only after dimensionless times + 
of the order of NV. To the accuracy with which we can now estimate NV, 
this may be a few decades. Thus, when a large change in ground water 
usage habits occurs (e.g. irrigation), all the the effects may not be evident 
for several years. Taking advantage of the foregoing, (4.2) is replaced by 
S, = (2/7)cAS,,, + (ey —vcos wt)S,,. (4.15) 
In the system in which the coordinates are tand v’ = v—(v/w)sin(wt —¢), 
this equation becomes 


II 


S, = (2 7)vAS, ¥ +ey'S,.. (4.16) 
Finally, for very large times, the solution becomes time independent, and 
S, in (4.16) can be replaced by zero to give 

(2/r)vAS,,,+ey'S, = 0. (+.17) 


5. A TWO-DIMENSIONAL MODEL FOR MIXING IN 
THE GUYBEN-HERTZBERG LENS 

The results of §+ imply that in a one-dimensional array of cells in which 
the fluid velocity distribution is ev+vcoswt, the mixing which occurs is 
equivalent to that which would occur in a fluid of diffusivity 2v\/7 which 
moved with the velocity ev relative to a coordinate system translating at 
speed «coswt. We adopt this result in order to simplify the analysis of the 
problem in which the salinity depends on the two space variables x and y. 
The actual Guyben-—Hertzberg lens has a mixing zone whose thickness 
decreases inland since the tidal response, and hence the effective diffusivity, 
decreases inland. We adopt coordinates such that y = 0 at the (translating) 
nominal interface position; this is essentially the y’ of §4. The distance 
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inland from the shoreline is —x. Let V(x)coswt—ye(x) be the vertical 
velocity distribution; let U(x) be the steady horizontal velocity; and let 
a(x) be the effective diffusivity 2)'(v)\ 7. In general V(x) will increase 
with increasing x, and conservation of mass requirements imply that 
U’(x)—e(x) = 0. With this notation the two-dimensional generalization 
of (4.17) is 


ta 
_ 
— 


U(x)S,.—e(x)yS,, = o(x)S,,,. ( 


The solution of this equation under the now familiar boundary conditions, 


S+lasy> «x, S+0O0as y > —&, is of the form 


S = }[1+erf{y/h(x)}], (5.2) 
and substitution of this into (3.1) shows that 
[h?(x)]' + (2e/ U)h? = 40 U. (5.3) 


For givene(x), U(x), o(x), and fora given h(x,) this equation can be integrated 
easily. 

Preliminary experiments on Maui island indicate that the physical 
facts are consistent with our predictions; however, there are too many 
guesses involved in choosing A, «, o, V to make any serious claims until 
more extensive measurements yield values for these quantities. 


6. THE EFFECTS OF MOLECULAR DIFFUSION 

Once fluid has been transported across the passage connecting two cells, 
the new fluid from cell »—1 and that already in cell » will be thoroughly 
mixed by the diffusion—mixing mechanism supplied by the irregular vortex 
motion in the cell. However, a few estimates are still needed to ascertain 
whether the transport of salt across the connecting passage is contributed 
primarily by convection or by molecular diffusion. ‘To answer this question 
we first note that the etfective diffusivity of the model of equation (4.1) is 
2V Az. (This is implied by the solutions which we have presented.) 

We now turn to the molecular diffusion process. We imagine the passages 
between cells to be holes of area a?, and assumea < \. The salinity gradient 
near and in the passage will be characterized by the quantity 6S a (8S is 
S,,—S,_,), and the total rate of flow of salt per unit time from cell n—-1 
into cell by diffusion (without convection) is characterized by vaédS, 
where v is the molecular diffusivity. Without a barrier, the salinity gradient 
is of order 5S A, and the diffusive transport of salt across the area \? is of 
order vSSA. The latter rate would be predicted by a conventional diffusion 
equation with diffusivity ». The former rate vbSa would be predicted by 
a conventional diffusion equation with diffusivity va’ \. We call this latter 
quantity the effective diffusivity, i.e. the diffusion coefficient which must 
be used in the diffusion equation to obtain the proper ditfusive transport 
for problems having the cellular geometry. The diffusivity which is 
associated with the mixing mechanism discussed in the preceding sectiéns 
of this paper is 27, \ + = VA. Thus, the relative effectiveness of these two 
transport mechanisms is given by the ratio va VA. 
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The value of a for a material of known permeability can be estimated 
in the following manner. The energy dissipation rate E per cell (assuming 
two holes to a cell) is of the order px( V A?/a*)?a*, where p is the fluid viscosity, 
VA? a? is the characteristic velocity near the holes and V’A?/a? is the velocity 
gradient near the holes. Furthermore, the definition of the permeability k 
states that gradp = ufV/k, where f is the porosity, and A’V gradp is 
the energy loss rate per cell since the pressure drop across the cell A grad p 
acts on a fluid area A? against a velocity V. Equating these two estimates 
of E, we obtain pV*Ata? = A8V2uf/k or kf} = a3/A. In the Hawaiian 
volcanic structure, kf~! is of the order 10-° cm?; and near the test wells 
on Maui, A can be estimated at 0-5 cm, v = 1cm?/day, V = 30cm/day. 
Using these figures, the effective molecular diffusivity va/A is approximately 
1/30 cm?/day, whereas 2VA/7 = 15. Thus, if the values of the parameters 
involved in any particular problem lie anywhere near the foregoing, the 
molecular diffusivity cannot play a competitive role in the determination 
of the salinity distribution. 


7. CONCLUSION 

In view of the lack of quantitative information regarding cell size and 
detailed velocity fields in permeable islands, we cannot quantitatively 
describe the salinity distribution in particular areas. ‘The observations 
that we do offer are that the discretized model of equation (2.2), the 
continuum model of (2.4) and the ‘intermediate’ model of (2.3) all imply 
the same quantitative predictions regarding the structure of the brackish 
zone for which the tidal fluctuations provide the mixing mechanism. The 
more highly simplified models introduced in §4 and §5 to deal with more 
complicated flow configuration are equally appropriate for the prediction 
of such salinity distributions. It is also clear that the predictions of the 
last of these can be evaluated without difficulty for any given porous medium 
properties and a given velocity field. 


Most of this work was carried out during a visit to the Scripps Institute 
of Oceanography in 1955. The author is indebted to W. H. Munk and 
Doak Cox for suggesting the investigation and for supplying the background 


informatton. 
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SUMMARY 
A mechanism is proposed by which cellular convective motion 
of the type observed by H. Bénard, which hitherto has been 
attributed to the action of buoyancy forces, can also be induced 
by surface tension forces. Thus when a thin layer of fluid is 
heated from below, the temperature gradient is such that small 
variations in the surface temperature lead to surface tractions 
which cause the fluid to flow and thereby tend to maintain the 
original temperature variations. A small disturbance analysis, 
analogous to that carried out by Rayleigh and others for unstable 
density gradients, leads to a dimensionless number B which 
expresses the ratio of surface tension to viscous forces, and which 
must attain a certain minimum critical value for instability to occur. 
The results obtained are then applied to the original cells described 
by Bénard, and to the case of drying paint films. It is concluded 
that surface tension forces are responsible for cellular motion in 
many such cases where the criteria given in terms of buoyancy 

forces would not allow of instability. 


INTRODUCTION 

Experiments have shown* that drying paint films often display steady 
cellular circulatory flow of the same type as that observed in the case of 
Huid layers heated from below. In the latter case (that of the so-called 
Bénard cellst) the motion can usually be ascribed to the instability of the 
density gradient that would be present if the fluid were stationary. This 
cannot be the mechanism causing the flow in the former case, since the 
circulation is observed whether the free surface is made the underside or the 
topside of the paint layer, that is, even if the gravity vector is effectively 
reversed. Instead it will be shown that surface tension forces are sufficient to 
cause instability and are probably responsible for many of the cellular patterns 
that have been observed in cooling fluid layers with at least one free surface. 


* These were brought to my notice by Dr R. Cousens of the Research Department 
of I.C.I. Ltd., Paints Division, who put forward the basic idea of a phenomenon 
dominated by surface tension, and pointed out many of the relevant physical factors 
that are introduced in this paper. The analysis given stems largely from his earlier 
investigation. 

+ For a description of these see, for example, Prandtl (1952). 


F.M. 2I 
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The essential physical features of the phenomenon spring from the 
fact that surface tension—in most fluids at most temperatures—is a monotone 
decreasing function of temperature and in the case of two constituents, 
a function of relative concentration. Thus if the free surface of a fluid 
layer is not at a uniform temperature or of uniform relative concentration, 
effective surface tractions are present and motion within the fluid must be 
expected to take place. ‘This idea is a very old one and has been used in 
a qualitative sense by several authors to explain many otherwise puzzling 
phenomena. 

If we consider for a moment the case of a homogeneous liquid layer 
cooling by radiation from the upper free surface and heated at the lower 
fixed surface, a simple qualitative explanation of the existence of steady 
cellular motion can be given. In the centre of the cells warm fluid is drawn 
towards the surface; this spreads across the surface, cooling as it does so 
until it reaches the edges of the cell where it descends towards the lower 
surface of the layer, and is there warmed. ‘The decrease in temperature 
across the surface from centre to edge of cell is accompanied by an increase 
in surface tension and hence by surface tractions that tend to maintain 
the circulation. The amplitude of the motion will of course be determined 
by the physical parameters that characterize the particular fluid in question 
and by the temperature gradients involved. ‘The driving force for the 
motion is provided by the flow of heat from the heated lower surface to the 
cooled upper surface. 

As in the case of Bénard cells induced by density variations, certain 
minimum requirements must be satisfied in order that cells may develop 
under the action of surface tension forces. These are examined analytically 
by means of small-disturbance theory (analogous to that developed by 
Rayleigh, Jeffreys and others for buoyancy forces) in the following section. 
Critical values for a certain dimensionless number are derived, corre- 
sponding to the case of marginal stability. Associated with these critical 
values are critical wave-numbers which define the size of the marginally- 
stable cells. 

In the concluding section an application of the results obtained to 
experimental observations on thin liquid films confirms the relevance of 
surface tension forces in many cases of instability. Somewhat surprisingly, 
it seems almost certain that many of the cells observed in molten spermaceti 
by Bénard himself were in fact caused by surface tension forces. 


SMALL DISTURBANCE ANALYSIS 

A small disturbance analysis will now be carried out for the particular 
case of an infinite homogeneous liquid layer of uniform thickness whose 
lower surface is in contact with a fixed plane and whose upper surface is 
free. The only physical quantities that are assumed to vary within the 
fluid are the temperature, the surface tension, which is regarded as a function 
of temperature only, and the rate of heat loss from the surface, also a function 
of temperature only. These are in fact the dominant factors in the cases 
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observed and a discussion of the relevance of this idealization will be given 
where required. For the case of a drying paint film, temperature variations 
can be replaced by concentration variations; the same analysis applies 
with only trivial redefinition of the parameters involved. 

The treatment adopted will be closely parallel to that developed by 
Rayleigh (1916), Jeffreys (1926), Low (1929), and Pellew & Southwell 
(1940), for the case of density-induced instability, and many of their 
arguments will be taken over directly. (A good concise account is given 
by Lin (1955, §7.3).) 

We start by defining the steady-state (unstable) equilibrium conditions 
that are to be perturbed, and we consider therefore a liquid layer whose 
bottom surface is given by the plane y = 0, whose free surface is given by 
the plane y = d, and whose temperature is a function of y only. If 

Oy, = rate of heat loss per unit area from the upper free surface in the 
unperturbed system, 
k = coefiicient of heat conduction in the liquid, and 
T\,,, = steady-state temperature of the lower surface in the unperturbed 
system, 
then the unperturbed temperature distribution 7), in the liquid will be given 
by 
Ty = Typ-— By, (1) 
where 
= kB, (2) 
since the rate of supply of heat to the surface from the liquid must equal 
the rate of loss of heat from the surface to its upper environment. In the 
relation (2), we suppose the magnitude of QO, to be defined by the upper 
surface temperature 7,. and the nature of the environment; the method 
of heat transfer may be conduction, convection or radiation, or any 
combination of the three. We shall assume that the environment exerts 
no mechanical forces on the liquid, and serves only as a heat sink. 

The temperature 7,, will itself depend on the nature of the liquid, 
the supply of heat to the bottom surface and the (equal) loss of heat from 
the top surface. However, we are not concerned here with a thorough 
analysis of the mechanism of heat transfer to and from the liquid layer, 
though these matters become relevant in the investigation of any particular 
physical phenomenon. Equations (1) and (2) provide sufficient information 
for our immediate purposes. 

Next, we superpose an infinitesimal disturbance and linearize the 
equations of motion and heat conduction. We write 

v = kinematic viscosity of the liquid, 

« = thermometric conductivity of the liquid, 

p = density of the liquid, 

t = time, 
velocity (supposed small) in the y-direction, 
= surface tension of the liquid, 


a 
II 


~H 
| 


212 
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O = rate of heat loss per unit area from the upper free surface (now a 
function of x and 2), 

T = temperature distribution in the liquid (a function of x, y and 2). 

T’ = T—T, is thus the perturbation temperature (also supposed small), 
and we write 7% for its value at the surface. 

The equations of motion and heat conduction become 


A 


(= " rv) 2p = 0, (3) 
_ cVi1¥'T’ = Bo. 4 
& 7 ) “) 


It will be noted that the equation for v contains no buoyancy term and is 
thus much simpler than in the case considered by Rayleigh. We also write* 

S = S,-—oT, (5) 
where —o = (eS/cT), Tog Tepresents the rate of change of surface 
tension with temperature, evaluated at temperature 7), with Sy = S(T7,,), 
and 

O= 0,+¢T;, (6) 
where g = (cO/cT)7 — 7, represents the rate of change with temperature 
of the rate of loss of heat per unit area from the upper surface to its upper 
environment. Qy is as defined before. o is a function of the liquid only, 
while gis likely to be affected in a complicated way by the surface-environment 
relations. ‘The relations (5) and (6) for S and QO are taken to depend linearly 
on 7” because we are considering an infinitesimal disturbance theory and 
therefore need only the first two terms in a Taylor expansion. 

The boundary conditions on the velocity are 


ov u 
o= = = ¥, (7) 
Cy 
when y = 0, 
es 
ey sonia 
z 0), pv=— = avi | igh (S) 
when yv = d, where 
4.) 02 2 
7 = A 
Cx" oz" 


Ihe second relation in (7) follows directly from the continuity condition. 
The necessity for the first of the relations (8), particularly in the case of 
neutral stability, is explained by Jetfreys. “The second relation (8) equates 
the change in surface tension due to temperature variations across the 


*'The boundary conditions (including (5) and (6)) are of crucial importance; 
by means of a suitable choice for these many physical phenomena may be very reason- 
ably idealized. The aim in this account is not to provide an exhaustive description 
of these phenomena and their relevant idealizations, but rather to provide a general 
treatment that illustrates the fundamental surface tension mechanism and compre- 
hends its many realizations. It is shown that the use of but three parameters suffices 


to describe the system. 
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surface to the shear stress experienced by the fluid at the surface, use being 
made of the continuity condition. The effect of the surface tension on the 
normal stress condition is ignored because its effects are small. (This same 
assumption is made in the conventional thermal convection problem. ) 

It is the boundary conditions on the temperature that need most careful 
explanation. It is worth considering in this respect the possible physical 
circumstances leading to the thermal equilibrium expressed by (1), (2), in 
order to understand the significance of the boundary conditions that will 
be applied to the perturbation temperature. ‘The supply of heat to the 
bottom surface of the liquid may be from a material whose coefficient of 
heat conduction is either large or small compared with that of the liquid. 
If the boundary consists of a good conductor of large extent, then the 
temperature 7, may be taken as a given fixed constant, and will lead, in 
the perturbed case, to the ‘conducting’ condition 7” = 0 at the lower 
boundary. If however, the bottom boundary of the liquid consists of a 
layer of a bad conductor itself in contact at its lower surface with a good 
conductor (for example, glass on metal), then the temperature gradient 
in the material bounding the liquid may be large compared with the 
temperature gradient in the liquid itself; 7), will, in the steady unperturbed 
state, be a constant, though it will be determined by the thermal properties 
both of the badly-conducting solid and of the liquid, and of their respective 
thicknesses. A slight change in the bottom surface temperature will not 
atfect the rate of heat conduction through the bad conductor by any 
appreciable amount. Hence the lower boundary condition on the 
perturbation temperature will be approximately c7’/cy = 0. ‘This we shall 
define, illogically perhaps, but consistently with previous authors, as the 
‘insulating’ boundary condition. ‘The two cases chosen above are limiting 


cases of the general mixed boundary condition 
T’=Y—, (9) 


which we can assume applies, where ¥ is a constant depending on the 
thermal properties of boundary and liquid. By a suitable choice of materials 
all values of Y ranging from 0 to « may be obtained. For convenience 
we shall consider in detail only the two extreme values Y = 0 and Y! = 0; 
since we shall find that they do not lead to radically different results, an 
exhaustive solution for arbitrary values of Y will not be attempted. 
Because of some confusion which has arisen among previous authors it 
is worth pointing out that the case ~' = 0 does not correspond only to 
a physically quasi-steady system in which the true insulating boundary 
condition, c7T/cy = 0, strictly applies, but also to a strictly steady system 
in which the so-called ‘insulating’ boundary condition on the perturbation 
temperature is a limiting approximation to the exact boundary condition. 

The balance between heat supply to and heat loss from the upper 
surface may exhibit similar characteristics. If, for simplicity, we consider 
a discrete jump in temperature as occurring at the free surface, then this 
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jump may be small or large compared with the drop in temperature across 
the liquid layer, depending on the efficiency of the processes for removing 
heat from the surface. Whatever the process, the equality 


~~ = gT" (10) 


must hold at y = d, using the relation (6) and the reasons given to justify 
the equality (2). If the temperature jump is large, then we expect q to be 
small, and vice versa. The particular value to be ascribed to q will depend 
on circumstances, and for this reason we shall retain it as a parameter in 
the subsequent analysis. It is worth noting that if g > oo, then we approach 
the ‘conducting’ boundary condition, and surface tension effects will not 
be present, whereas if g->0, we approach the ‘insulating’ boundary 
condition, as defined above, valid for perturbation of a strictly steady-state 
system. 

A further application of this small-disturbance analysis to concentration 
gradients in a layer composed of a mixture of two liquids will show—if 
only by comparison with experiment—that predictions of instability are 
relevant for quasi-steady systems. If convection cells can be set up in 
a time small compared with that in which changes occur in the quasi-steady 
thermal gradient §, then the application of strict ‘insulating’ boundary 
conditions leads to meaningful results. 

If we now introduce dimensionless variables 








x 3 tk 
E,n,¢ NG rs ie iia. 
ee? (5 d 5) d aad 
and suppose that 
v= —<F(E,o)f(n)e! (12) 
d 
T’ = Bd F(E, C)g(n)e”’, (13) 
we find that 
— = +o2F = 0, (14) 
Of as if Say 


where x is a dimensionless constant arising from the separation of variables ; 
and that 


[ p — Pr(D? — x?)](D? — x2) f = 0, (15) 
[p — (D?—22)]g = -f, (16) 


where D = d/dy and Pr = v/«x. The particular choice of the forms (12), (13) 
for v, T’ is well explained by Pellew & Southwell (1940) who show that 
the constant x is merely indicative of a periodic structure in the (€, ¢)-plane. 
The shape of the cells associated with the solution obtained is not specified 
and a second-order theory would be required to select a particular cellular 
structure. Christopherson (1940) has given the solution of (14) for 
hexagonal cells. 
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The boundary conditions (7), (8) and (10) become* 


FO) =f'(0) = 9; fl) =9, f°) = #Ba(1), g’(1) = —Lg(1); (17) 


where 2 
pu SE, (18) 
pvkK 
and qd 
L== 9 
: (19) 


are dimensionless constants. Physically speaking, the number B is a 
dimensionless grouping of physical parameters expressing the relative 
importance of surface tension forces (caused by temperature variation) 
and of viscous forces, in any small disturbance responsible for both of 
them. The number L can be suitably interpreted in terms of the arguments 
given earlier with regard to relation (10). 

As a last boundary condition, we have for the conducting case (Y = 0 


in (9)), 9(0) = 0, (20) 
and for the insulating case (Y~! = 0 in (9)), 
g'(0) =0. (21) 
The general boundary condition (9) would lead to the relation 
g'(0) = Mg(0), 


where .W = Y/d, and although we shall not retain it as an independent 
parameter, we see that ./ forms with B and L a set of three parameters 
which suffice to describe the system for our purposesf. 

In common with earlier workers, we look for solutions of the case of 
neutral stability when we put p= 0. Equations (15) and (16) become 


(D? — z?)(D? — x?) f = 0, (22) 
(D?— 22)g = f. (23) 
The solution of (22), subject to the first three conditions (17), is 


x cosh x — sinh x 





f = a( sinh xy + =e 7 sinh xn — xy cosh a) (24) 


and the solution of (23), subject to the last condition (17) and (20) or (21) 





becomes 
a xcoshx—sinhz , 

g = al —7coshzy + 7? cosh an — 
4x, 4x sinh x 


xcosha—sinhx . 
7 sinh xy — 





— 47? sinh xy — 
4? sinh x 


(x? cosh? x + % sinhzcosh x + sinh? x) + L(«? + «sinh z cosh + sinh?«))} 
ae a x 


42? sinh x(x cosh x +L sinh z) 





as 


x sinh =| (25) 


*(’) now denotes differentiation with respect to the single independent variable. 

+ It must be made clear that the evaluation of these three parameters in any physi- 
cally observed circumstances is not necessarily easy; it is however a separate problem. 
A lengthy consideration of possible values to be expected, beyond that already given, 
would only obscure the explanation of the fundamental instability that is being 


analysed. 
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for the conducting case, and 


3 xcoshzx—sinhz , 
g = a| — 7cosh27-+ 
7a 4 sinh x 





x.cosh x — sinh x 3 








4x? sinh x 


(xcosh x — sinh z)* + L(x? + xsinh x cosh x —2 sinh? x) : 





4x? sinh x(x sinh x + Z cosh x) 

(26) 

for the insulating case. Substitution into the fourth of the boundary 
conditions (17), namely the one involving both f and g, vields a relation 


between B, L and x. This ts 


Sx(z% cosh x+Asinh z)(z— sinh xcoshz) “ 





(x? cosh x — sinh? x) 


| 


1 . : 
for the conducting case, and 


R Sx(xsinhz+hcoshx)(x—sinh xcoshz) a 
) ens eeaememe. ars — (25) 
, awe nn 
(x® sinh x — x* cosh xz + 2zsinhz-— sin! 








tor the insulating case. These are plotted in figures 1 and 2 for the various 











values of L. The curves obtained are neutral stabilitv curves, and without 
repeating the analysis tor p # 0 it is clear which regions refer to growing 
| | i ' ‘<3 eo et 
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ey ] 1 , ae 17 ee noniivn ] ’ + 
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condition is particular in that the critical value of z is zero, B = 48. Apart 
from this special case, the shape of the neutral stability curves 1s very similar 
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would depend ¢ e thickness of the nit and I verv th hims ould 
tend to zero Det ed solutions for the more gene ( rv condition (Y) 
— nor rer i. 7 | re rt + wane 9 bass alee ae ner } 
with finite non-zero values of } are not given Dut it ls not expected that 
‘. 1d 1 ] - iokwean ae } 7 ew . 
these would behave verv differently from the two limiting cases } (), 


t mnare the dimencant number B. defined I 1&8 
tO COMpare the GlmMensiloniess numoe >, denned DV (19), 


hich : r r he srt: ) > y thanicr s1-31th ) ) ol .. 
which 1s relevant for the surface tension mechanism, with the Ravleigh 








1 1 ] j } roa = y ¢} asin 
where g is the gravitational acceleration and » is the coefficient of thermal 


; 3 ] 2a which 3 le st tor th} lane} j nde ,eECh 
expansion of the liquid, which ts relevant for the density-dependent mech- 


anism. In particular we may compare the critical value B.., = 80 when 
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L = 0, from figure 1, with the corresponding critical value 4. 


. = 571* (see 
Jetfreys (1926), the case of one free ‘insulating ’t and one fixed conducting 
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surface, and the amending paper Jeffreys (1928)). The former of these 
leads, for any given fluid layer, to a critical thickness 


80 pvne\12 
d, =( | 
op 


and the latter to a thickness 








These will be equal for a value d,, given by 


ad, = — —. 
“ — pgy 80 

For most liquids at laboratory temperatures this relation leads to values 
for d,, of the order of lcm. For thicknesses less than d),, then, we expect 
surface tension forces to be more effective than buoyancy forces in producing 
instability, and, for values of d as small as 1 mm, the onset of cellular 
motion could confidently be attributed to surface tension rather than to 
buoyancy. Since L will not be large for sufficiently small values of d and 
since B.,,, does not vary very rapidly with L, the arguments are not altered 
greatly if non-zero values for LZ are assumed to apply. 


COMPARISON WITH OBSERVATION AND EXPERIMENT 

The original observations of Bénard (1900, 1901) provide a very suitable 
illustration of the instabilities considered above. For with spermaceti 
he achieved cells in layers less than $} mm thick when the temperature 
of the lower surface was maintained at 100°C. Moreover, these cells 
persisted even when the temperature of the lower surface approached 50° C. 
Since spermaceti solidifies at rather above 50° C an upper value can be placed 
on the temperature variation across the layer and can only have beena degree 
or two, if that, in some cases. Bénard gives a value of 7-7 x 10~* c.g.s. units 
for y and, fortunately for our purposes, a value of nearly 5 x 10~* c.g.s. units 
foro. We can assume that the case of a lower rigid conducting and an upper 
free ‘insulating’ boundary applies approximately to his experiments since 
the lower surface was metallic and the temperature difference between the 
upper surface and its surroundings can be assumed large compared with 
the drop in temperature across the liquid layer. It would require detailed 
calculation to obtain an accurate estimate for L in this case, since the 
mechanism for heat transfer involves both convection and _ radiation. 
However, a rough calculation shows that L < 1—whether it be 10-3 or 107! 
is scarcely important since figure 1 shows that the critical value B necessary 
for instability only increases by a factor of 3 when L = 5—and we can 
therefore use the approximation L = 0 as an initial approximation. 

Even if we assume that « was as small as 10-4 and v as small as 10, 
the value obtained for 4 is still less than 571. On the other hand, for 
reasonable values of « (5 x 10-4) and v (5x 10-?), a value of B in excess 
of 80 is obtained even for a temperature drop of only 1° C across the layer. 
Furthermore the ratio of cell size to layer thickness quoted by Bénard 





On convection cells induced by surface tension 499 


leads to values of « lying between 2-1 and 2-5 which are in good agreement 
with the analysis given in the previous section, where «,,,, = 2-0. The 
equivalent value of «,,,, from the analysis of Jeffreys is 3-5, which even 
allowing for error is significantly different from the observed values. 
Thus we see that the buoyancy mechanism has no chance of causing 
convection cells, while the surface tension mechanism is almost certain 
to do so, and that observations support this. (Going back to the 
choice of L = 0, we see that this was not critical and hence that an exact 
analysis of the heat transfer at the surface is not necessary to sustain the 
argument.) An intimation that the instability theory based on buoyancy 
forces would not account for all of Bénard’s results appears in a paper by 
Volkovisky (1939). 

The allied problem of a liquid cooling by evaporation may be treated 
in a similar fashion. Although no strictly steady equilibrium state may 
be presupposed, owing to the loss of fluid from the surface, the effect of 
evaporation may be reasonably well represented by a given heat loss from 
the surface, using a value for L that depends on the rate of evaporation 
and the latent heat of vaporization. A quasi-steady value for 8 must be 
assumed to hold. 

The application of the analysis to the case of a mixture of two liquids, 
one more volatile than the other, requires more explanation. ‘Iwo factors 
tending to instability will now be relevant. ‘The first, that due to temperature 
variations, has been treated already. ‘The second is that due to relative 
concentration variations. ‘This may be treated separately using a suitable 
interpretation of the analysis given above for temperature effects. For 
if the temperature TJ is interpreted as the concentration C of the volatile 
liquid /, in the non-volatile liquid /,, the constant & is interpreted as the 
coefficient of diffusion of #, in #, and the parameters o, g and y refer to 
variations with concentration C, then the equations (3) and (4) are 
appropriate equations for describing small disturbances due to concentration 
variations in a layer whose equilibrium (or quasi-steady) concentration 
gradient is given by B—for the moment we regard 8 as some unspecified 
function of y. Previous work has shown—see for example Batchelor (1954) 

-that the most unstable case, for any given concentration drop across the 
layer, occurs when the gradient f is a constant*. Consequently, although 
the instantaneous concentration across the layer for the casé of surface 
evaporation of /, from a mixture of #, and #, would never be a linear 
function of y, we may safely use a value of 

B= C.../d 
to obtain a lower limit for stability. The boundary conditions on the 
velocity, (7) and (8), would still apply; the relation (10) would also apply 
with suitable interpretation; but at the bottom surface the insulating 
condition would now be relevant. Hence the result (28), displayed in figure 2, 
can fairly reasonably be applied to describe the concentration effect. In 


* This has been verified analytically for the particular choice of B: B = C,/bd, 
d>y>di-—b);8 =0,d1—6b) >» >0;0 <6 <1. Instability is greatest for b = 1. 
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general terms the temperature and concentration effects are additive and 
for any given circumstances a joint criterion replacing either (27) or (28) 
could be calculated. ‘This will not be carried out here because of the 
difficulty in prescribing the equilibrium concentration gradient. 

A series of experiments carried out in this laboratory on thin layers 
of mixtures of liquid paraffin and ether have shown convection cells of 
polygonal type when evaporation of the ether takes place. Estimates of 
the buoyancy forces involved showed that these could not be responsible 
for the motion. On the other hand, the difference in surface tension 
between liquid paraffin (32 c.g.s. units) and ether (19 c.g.s. units), combined 
with their temperature variations, leads to surface tractions large enough 


0 


to maintain the motion. For example, a 5°, concentration of ether 
evaporating at room temperature (say, 15°C or less) showed cells down 
to a layer thickness of 0-5 mm, this critical thickness being highly repro- 
ducible; even allowing for a very high rate of evaporation, and adding 
the concentration and temperature effects on the density, a value of 4 
of the order of 10 or less is obtained. Similarly crude estimates of B for 
the two effects suggest instability. ‘The ratio of cell size to layer thickness 
proved to be relatively constant over a fair range of layer thickness and 
using the analysis of Christopherson (1940) a value of z ~ 2-0 was obtained, 
which is in reasonable agreement with our analysis. 

Similar experiments using a wide range of both volatile and non-volatile 
liquids have been carried out by Dr Cousens. All his results are in excellent 
qualitative agreement with the ideas described here. In particular, by 
choosing a volatile fluid with higher surface tension than the non-volatile 
base, cellular motion was inhibited. ‘The addition of a very small quantity 
of silicone oil to various mixtures also inhibited cell formation by reducing 
the surface tension overall. 


‘The author wishes to thank Dr R. Cousens for pointing out the probable 
role of surface tension and for helpful discussion of the observations quoted. 
It is hoped to publish a joint paper on the specific application of the 
foregoing theory to paint films. The author ts also grateful for many useful 
criticisms by the referees. 
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A note on shock flow in a channel 


By ROY GUNDERSEN 
Department of Mathematics, Illinois Institute of Technology 


(Received 20 December 1957 and in revised form 25 March 1958) 


SUMMARY 
A shock wave is generated by a uniform compressive piston 
motion and passes into a channel of slowly varying cross-section. 
A relation in closed form is obtained between shock strength 
and the area of the channel and is used to discuss converging 
cylindrical and spherical shocks. 


1. INTRODUCTION 

Using a linearized theory based on small area variations, Chester (1953, 
1954) discussed the motion of a uniform shock wave passing through a 
two-dimensional channel composed of two uniform cross-sections separated 
by a section of slowly varying cross-section. ‘The fluid in front of the shock 
was at rest; initially, the flow behind was isentropic, but when the shock 
entered the transition section the shock strength was altered and the 
subsequent flow was not isentropic. The disturbance depended on the 
area but not on the shape of the cross-section; it consisted of two 
perturbations: a permanent disturbance due directly to the area change, 
and a transient disturbance, propagated with sonic velocity relative to the 
flow behind the shock, due to the reflection of the permanent disturbance 
from the shock. The pressure change behind the shock was determined. 
Chisnell (1957) integrated this first-order relation to obtain a functional 
relation between the channel area and the shock strength, and used particular 
channel shapes to discuss converging cylindrical and spherical shocks. 
His results agree closely with previous similarity solutions obtained by 
Guderly (1942) and Butler (1945). 

In Chester’s solution, the shock has, so to speak, come from infinity, 
and there is no possibility for reflections to occur upstream of the shock. 
In a discussion of the non-steady flow in a channel of slowly varying 
cross-section, linearized with respect to small area variations, the present 
author obtained Chester’s solution on the basis of one-dimensional theory. 
A simple but interesting generalization is to assume that a shock is produced 
by a uniform compressive piston motion; when it passes into a channel of 
varying cross-section, there are three distinct contributions to the 
disturbance: a permanent disturbance due to the area change, a transient 
disturbance due to reflection from the shock, and also a transient disturbance 
due to reflection from the piston. ‘This problem was solved by the author 
(Gundersen 1958) for the case of a slowly converging or diverging 
cross-section, and an expression was given for the pressure change 
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directly behind the shock. ‘This differential relation is integrated here 
in closed form to give a functional relation between channel area and shock 
strength. Since the shock strength is uniform over its area for the 
symmetrical converging cylindrical and spherical flows, this result can be 


utilized to discuss such flows. 


2. "THE RELATION BETWEEN CHANNEL AREA AND SHOCK STRENGTH 
The fluid is assumed perfect with constant specific heats. Denote 
by y the adiabatic index, by Py, and c, the pressure and sound velocity 
in the gas at rest in front of the shock, by w the shock velocity, and by ws, 
P, and c, the fluid velocity, pressure and sound velocity behind the shock. 
Let w, be the shock velocity relative to the flow behind it, so that w, = w— Us, 


and let M,, 7, and m denote the following Mach numbers, VM, = wc, 








M, = @y/€s, m = Uy/cy. ‘Then from the usual formulae for a normal shock : 
u y-1+2M;? 
Ws 1+2M-? M? = — —— , 
a= ( eres 2y—(y-1)M* 
2w(1— M7?) 2(1—M?) — P. 
2 a, = ’ »» = (2yM?2-y+1 : 
“s y+i m= G+ lM; 2 = ( eres 


It is convenient to express the cross-section of the tube in the form 
E(x) = E,+E,(x), where E, is the original uniform cross-sectional area. 
In terms of M,, M, and m, the pressure perturbation directly behind the 
shock, for E,(x) = kx (where & is a constant, positive for a diverging section 
and negative for a converging section) and when the shock is produced by 
a uniform compressive piston motion, is 

P, = —(P,—P,)E, EK. (1) 
There is a small error in the previously presented A, which is due to the 
omission of a factor of 2 in the next to last equation on page 564 of the 
previous paper, which should have read 


Bo _ of uo— (1-0), 
Po af W (tty — W,) «( ‘ 


so that equation (8.2) should have read 
. (Ko [uy—(1-4#)Wo] 
Sp = 2c fy—1)4— —| a | Se( 2). 
io = Ded IM | ay |e 
This last expression also appears on page 575 with different subscripts, 
and the correct expression in terms of the parameters M,, M, and m is 
$2 = Cyo(y—1)w-h{4(1 — Mi) [(y +1) - My) 8+ 
+m? —2+2(y—1)(y—1+2M>?)—}}. 
When this expression is substituted into the defining equation for K, on 
page 576, the correct K of (1) is defined by 
2(y-—14+2M>*)(y+1)7K- 
= 2M?2+mM,(1— M;?)+(1+M,;*)—(y—1+2M>5*)(2y)-? x 
x {4(1 — M7) [(y + 1)(1 — Mp?) 1 +m? - 24+ 
+2(y—1)(y—1+2M5*)}, 
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and is a monotonic decreasing function of the shock strength. For y = 1-4, 
K varies from 0-500 for weak shocks to an asymptotic limit of 0-259 for 
strong shocks. In Chester’s solution, the variation of the corresponding K 


is from 0-5 to 0-394. 


‘T’o facilitate comparison, Chisnell’s notation will now be utilized when 
possible. If the initial shock strength is defined by P, P,; = s, (1) takes 


the form: 











dé | 1] 

Eds (z—1)K(s) 
Ff - iie-T) (y—1)[(y+1)s+y-—1] 
2y?s[(y—1)s 1] 2y*s(z—1) 





a . (2) 
-1)[(y-1)s+y4+1] 


Integration gives the relation between area and shock strength, but as 
Chisnell points out, for the case of two uniform channels connected by a 


transition area of varying cross-section, the result gives the average strength 
of the shock after the area change only if the total change is small. For 
large area changes, it will be only an approximation to the average strength. 


The integration of (2) gives 


E f(z) = const., 


where 


fla) = (= 1)s0-% [y= Le ty + 1], 


For weak shocks (sz nearly unity), f(z) behaves like (s—1)? for y = 1-4, 
i.e. the acoustic result that the strength of the disturbance varies inversely 
as the square root of the area of the disturbance. 

As pointed out by Chisnell, the above results are applicable to converging 
cylindrical and spherical shocks if channel areas are proportional to R or R?, 
respectively, and R is the distance of the shock from its axis of symmetry. 

Near the point of collapse of the symmetrical shocks, z is very large, and 
f(z) behaves like z1*, where K, the asymptotic limit of A(z), is 


, 2y(y—1) 


(2y—1)(y+1)° 


Hence, for cylindrical and spherical shocks near their axes of symmetry, 
the strength is proportional to R-* and R~*4, respectively. The following 
table compares the values of K with those obtained by Chisnell for the 


corresponding number in his paper. 


Chisnell 
y= 0-326223 
y= Z 0-394141 
y=3 0-450850 


This paper 
0-155844 
0-259259 
0-357143 
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For y = 2, 4, 3, the numbers of the present paper are about 52%, 34% 
and 21°, respectively, smaller than the numbers given by Chisnell, i.e. 
the perturbations reflected from the piston reduce the ultimate shock 
strength, near the point of collapse, to a value less than in the Chisnell 


determination. 
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On the melting of a semi-infinite body of ice 
placed in a hot stream of air 


By LEONARD ROBERTS* 
Department of Mathematics, Massachusetts Institute of Technology 


(Received 7 February 1958) 


SUMMARY 

A simple mathematical model is proposed to describe the 
steady melting of a body of ice which presents a plane surface 
transverse to a stream of hot air; the temperature of the air is 
such that vaporization does not occur. 

The analysis takes into account the convection of heat away 
from the surface by the water released in melting and the results 
show that the rate of transfer of heat to the body and thus the 
rate of melting, is reduced by as much as 46°, by this convection. 

Simple approximate expressions are obtained for the rate of 
melting, the thickness of the water layer, and the thickness of the 
thermal boundary layer in the ice, in terms of a basic parameter S 
which can be calculated in terms of known quantities. These 
results are compared with those obtained by a separate Pohlhausen 
calculation and are found to be in good agreement 

It is also shown that there exists a thermal boundary layer, 
in the body, of thickness much greater than that of the boundary 
layer in the air, in which the temperature changes rapidly from 
its value at the melting surface to its value in the far interior. 


1. INTRODUCTION 

When a non-insulated body is placed in a stream of hot air a transfer 
of heat to the body takes place; when the temperature of the stream is 
high enough the body may melt and even vaporize. The rate at which 
the body melts is determined by the amount of heat available for latent 
heat; this is the difference between the total amount of heat transferred 
to the surface and the amount conducted away from the surface to the 
interior of the body. 

The present paper considers the steady melting of a semi-infinite 
body of ice whose plane surface is placed normal to a stream of hot air 
(see figure 1); the analysis may be regarded as applicable to the conditions 
in the neighbourhood of the forward stagnation point of a blunt-nosed 
body of arbitrary shape. ‘The temperature range considered is such that 
vaporization does not occur. It is possible, under certain simplifying 
assumptions, to reduce the Navier-Stokes and energy equations to two 
ordinary differential equations which are solved by an approximate method 
for both the plane and axisymmetric flows. 


* Now at Langley Aeronautical Laboratory, N.A.C.A. 
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An interesting and important feature of the flow is the shielding effect } 

of the thin layer of water which forms between the ice and the air. The 
rate of heat transfer to the body, and thus the rate of melting, is reduced 
in two ways by the presence of the water layer: firstly, the temperature 
of the air—water interface is raised above that of the melting surface, 
which tends to reduce the rate of heat transfer from the air, and secondly, 
convection of heat takes place in the water laver which further reduces 
the rate of heat transfer across the melting surface. 
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It is found that the first shielding etfect is small, but that the second 
shielding effect may become appreciable. In figure 2, there is shown the 


shielding ratio R,, defined as 


R rate of transfer of heat from water to 1ce 
x ne Ss ess ae SER? 


rate of transfer of heat from air to water 








wi 
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vs the basic parameter S defined by 
_ p(T. -T,,) 
L+ce,(T,,, —T_,) 

where c,,, is the specific heat of air, c, is the specific heat of ice, L is the 
heat of fusion of ice, and T,, T,,, and 7_, are respectively the temperature 
of air at large distances from the body, the temperature of water at the 
melting surface (taken always to be 0 C), and the temperature in the far 
interior of the body. (The temperature 7* at the air—water interface also 
plays an important role in the theory, but has been eliminated from the 


S 





final results.) 

It should be noticed that the shielding effect is largely independent of 
the flow conditions, since for ordinary conditions, the shear acting on the 
water layer is always sufficiently large to carry the water away at a reasonable 
speed compared with the rate of melting. On the other hand, the details 
of the flow in the water layer, such as its thickness and the rate at which the 
ice melts, do depend on the conditions of flow. 


2. DESCRIPTION OF THE FLOW AND METHOD OF SOLUTION 


2.1. The flow configuration 


The flow considered is that shown in figure 1. In order to reduce 
the problem to a steady state the coordinate system chosen has axes fixed in 
the melting surface; in this coordinate system the interior of the ice moves 
towards the stationary melting surface, s = 0, with constant velocity, w,,,, 
equal to the rate of melting. 

Considerations of the continuity of mass lead us to expect (under the 
assumptions of §2.2) a water layer of constant thickness between the ice 
and the air stream. From the Hiemenz solution for flow near a stagnation 
point (which will be recalled in detail later) is is known that the horizontal 
component of velocity in the air is proportional to the horizontal distance x 
from the stagnation point, and it follows that this is also true of the motion 
of the water since the velocity is continuous at the interface. Also, since 
the rate of heat transfer from the air, in the vertical direction, is independent 
of x, it is to be expected that the body will melt at a rate which is independent 
of x. Thus the amount of water crossing a vertical line at the position 
x = x, is that which has been produced by melting between the line of 
symmetry x = 0 and the line v = x,; this amount is proportional to x. 
The ratio of the amount of water to the velocity is constant, i.e. the thickness 
of the water layer is constant. 

Because of the continuity of velocity and stress at the air—water interface 
=< = s*, the motion of the air is affected only slightly by the presence of the 
water. (It will be shown that the velocity of the air at the interface is of 
the order of the square root of the density ratio, (p,,./p,.,,..)!2 which is 
O(10-*) for most of the temperature range under consideration.) The 
Hiemenz solution, modified to take account of the variation with 


2K2 
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temperature of the density, viscosity, and thermal diffusivity of air, is a 
good approximation for the flow now under consideration. In a region 
near the stagnation point the flow velocities are small and compressibility 
effects are ignored. 

In the water it is assumed that the density and the thermal diffusivity 
are constant but the viscosity varies with temperature. 


2.2. The basic assumptions 

Only the case of laminar flow is considered in this paper. In addition, 
the follewing approximations are made in order to simplify the analysis. 

(1) The magnitude of the components of velocity in the air are such 
that the approximation WM = 0 (M = Mach number) may be made. The 
equation of state is then p,7’ = constant, where p, is the density and T is 
the absolute temperature. 

(2) The viscosity ., and thermal conductivity k, of the air obey the laws 


~~ 


= constant, = constant. 


“1 
7; 7 





(3) The specitic heat of air c,,, 1s constant. 

(4) ‘The Prandtl number for air, o,, is constant (this is implied by (2) 
and (3)). 

(5) Viscous dissipation is negligible. 

(6) The density p,, thermal conductivity k,, and the specific heat c,, 
of water are constant, but the viscosity uz, varies with temperature. 

(7) The melting temperature 7,,, the specific heat c,, the thermal 
conductivity k,, and the density p,, of ice are all constant. 


m? 


2.3. The method of solution 

‘The solution is developed in the three regions, s > s* (air), = < 0 (ice) 
and < = < s* (water). Atthe interfaces s = s* ands = (0 certain quantities 
which appear in the boundary conditions are not known, a@ priori, but are 
part of the answer; this introduces a number of parameters which must 
be determined as part of the overall solution. The most important of these 
ire 7*, the air—-water interface temperature, 2*, the water laver thickness, 
ind w,,, the melting rate (or the non-dimensional counterparts of these 
quantities ). 

In the actual computation, families of solutions were found corresponding 
to a range of these parameters and the appropriate member of the family 
determined by matching the solutions at the interfaces. 

The only truly arbitrary parameters in the problem are T,, T_, and §,, 
which characterize the outside potential flow of the air (u ~ 8,s ass— x), 
and the solution must be expressed in terms of these three quantities. 
In some sections of the paper it is easier and more convenient to express 
the solution in terms of 7* rather than 7, but all results are given 


eventually in terms of 7,, T_, and §,. 
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The solution of the equation for heat transfer in the melting ice is first ~_ 


obtained in simple exact form, in terms of the unknown melting rate ~w,,. 
Attention is then turned to the equations which govern the motion of the 
air. ‘These are reduced to a system of two ordinary differential equations 
of which an approximate solution is obtained involving two unknown 
quantities, 7* and ¢ (a small parameter occurring in the velocity at the 
air—water interface). 

The equations for the water layer are also reduced to ordinary differential 
equations and it is shown that the number of boundary conditions to be 
satisfied at s = 0 and : = s* are just sufficient to determine the solution 
in the water layer and all the unknown parameters that are introduced. 

Simple estimates are then deduced, without solving the equations for 
motion in the water layer, which give useful approximate results concerning 
the rate of melting and the shielding effect of the water. Finally, a better 
approximation which uses the Pohlhausen method to determine the water 
layer thickness and the approximate velocity and temperature profiles in 
the water is made; this also gives a check on the first simple estimates. 


3. HEAT TRANSFER IN THE ICE 


The transfer of heat in a solid moving with constant velocity w,, (unknown) 


m 


in the z-direction is governed by 











C f C ( C T C C fi 4 
PzlgW,, —— = —| kg —])+ —( kg —}; (3.1) 
cz COX\ ox cz os 
the boundary conditions for the present problem are 
T =T at = = 0, | 
; (3.2) 


T=T as [> — ©, 
where the quantities ps, cs, k,; and 7',, are those defined in § 2.2. 
The solution of (3.1) which is a function of z only and satisfies (3.2) 
is obtained by elementary methods and is written as 


ww 


T=T,+(T,-T Jexp( i aa 
\ 3 


This gives the rate of heat transfer at zs = 0 to the interior of the ice, 


C d ia aa | 2 
(x: =) = tt tT. ~T.); (3.4) 


+ 


equation (3.4) will be required later in the formulation of the boundary 
conditions, at z = 0, for the equations which describe the motion of the 


water. 

The form of the solution (3.3) shows that there is a thermal boundary 
layer near the surface s = 0 in which the temperature changes from 7, 
to T_,,; the thickness of this thermal laver is O(k,/w,, p3¢3) and thus varies 
inversely with the rate of melting and directly with the thermal diffusivity, 
Rs /'p3 C3. 
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4. THE MOTION OF THE AIR 

4.1. The equations and the boundary conditions 

Consider the steady flow of air of variable density p,, viscosity 4,, and 
thermal conductivity k,, at zero Mach number in two dimensions. If the 
components of velocity in the x- and 3-directions are u and w respectively, 
and the pressure is p, then the equations of continuity, momentum, and 
energy are 


(piu) , are) 

















— =) (4.1) 
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u of + ey M4 ! ae pe (1 =) gp’ (1 atl ; (4.4) 
CX oz p,a, | ex OX oz\' “ oz 
Where c,, = constant and o, = 4,¢,, Rk, = constant are respectively the 
specific heat at constant pressure and the Prandtl number, and dissipation 
has been neglected*. 
\t large distances from the body the flow has the potential form 
“= By x, w 2, De Fins 


hat the boundary conditions for the problem under consideration are: 


t 
so t 


. . ae eRe ye 
i= el u 0), 1 1 ‘ at 2 Pe (4.5) 
un fix, T>T,, as 3—> OO. 
Here 8, and 7’, are given by the potential flow; 7* and ¢ (a small parameter 
which represents the etfect of the water layer on the motion of the air) 
are unknown and will be determined by the rate at which the ice melts. 

It is well-known that the Navier-Stokes equations for incompressible 
How have an exact solution which represents flow near a stagnation 
point. This solution, the Hiemenz solution (Goldstein 1938, vol. 1, 
pp. 139-140), has the form 

u = B,xf'(n), w= —(8,u,/p,)'?f(n), 

where 


“» 


ee (8, py 4)! 


“The terms y/pc,(@u ¢z)? and w/pc, @p/@x which represent the change in 
temperature due respectively to the heat generated by skin friction and the adiabatic 
compression are O((y—1).W*) compared with other terms in equation (4.4), and 
have been neglected. A comparison of the results of Levy & Seban (1953) and 
Brown & Donoghe (1951) shows that this neglect is justified. 








The melting of ice in a hot stream of air 511 


and f satisfies the differential equation 
f?— ff" =1+f" (4.6) 
and the boundary conditions 
f(0)=f'0)=0, f'(x)=1 (4.7) 
(a dash denotes differentiation with respect to 7). 

In this solution the fluid is assumed to have constant properties so that 
the temperature, 7 = 7)+(7,.—T ))g (7) = wall temperature) may be 
calculated separately from the reduced energy equation 

g’ +o,fe’ = 0, 
with 2(0) = 0, g( co) = I. 
For such a flow the boundary layer has constant thickness and the vertical 
component of velocity and the temperature are constant in layers at a 
uniform distance from the wall. 


(4.8) 


In an analogous way we seek a solution which has similar properties 
but which takes account of the effect of the temperature field on the density 
and viscosity (through the assumptions of §2.2) and consequently on the 
velocity field. 

We consider a solution of the form 





” pP of 
“= B, xf, (71), ee SS = (B; Px/P )" *fi(m), 
? Px (4.9) 
P Mey. . oo a | 
PF = 8iln)s r dy, = (Bi p/p)? (2 -2*), 
where p, T= p,,7,, and p,/T = u,/T,. With new dependent variables 


fi(m,), £:(7,) and the new independent variable »,, (4.1) is satisfied, 
(4.2) becomes 
fy -fi Si = Fi fy , (4.10) 
ind (4.4) becomes 
git+o fig, = 9, (4.11) 
where a dash now denotes differentiation with respect to 7). 
The pressure p, is given by 


p = —3p, Bi x*?—7(7,)8, b-, 
where, from (4.3), 7(7,) satisfies the reduced momentum equation 
mW =tAlh 2) +18h 2) +3 fil -afi ° (4.12) 
The boundary conditions (4.5) reduce to: 
T* 


f,(0) = 9, fi) =«, g,(0) = T.’ fi(2%) = 1, g(o)=1. (4.13) 


An approximate solution of (4.10) and (4.11) with the boundary 
conditions (4.13) is obtained by considering a perturbation, in ¢, of the 
particular solution for which e = 0, i.e. 


fi = $9 +€$,+ Ofe?), and g, = dy +ep, + O(e?). (4.14) 
The solution (do, %) is the first approximation and takes no account of the 
effect of the motion of the water. 
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4.2. The initial approximation, « = 0 

When « 0, the boundary conditions (4.13) are those for flow past a 
solid boundary; the approximate solution of the system (4.10), (4.11) 
with (4.13) has been given by Levy & Seban (1953); it was assumed that 


the solution ¢), y%) has the approximate form 


y/ ae Cs d 1 ) 
Dp = | - exp ay), + 77 ny + ra + 7 ut} 


; P ; (4.15) 
; Bones ; ) 
bs, 1 - Kexp( 4», + 7 "Ty + 31 1) + tI it). 


so that the boundary conditions as 7,-- & are satisfied if d and D are 
negative; the non-dimensional temperature at 7, = 0 is Y%(0) = (1—K). 
The constants a, 6, c, d, and A, B, C, D, are determined by satisfying the 
boundary conditions and the differential equations as far as possible at 
7, = 9 ‘The system of algebraic equations which results reduces to the 
following pair which give a and A in terms of K 
10K A[2a? —(1—K)] = 24a° —5a®(12(1—K)—1]-4 
+2a[15(1—K)*?-(1—A)], (4.16) 
244*®°—50,aA? = —o0,(1—K)A; 
the remaining unknown coefficients are then easily expressed in terms of K. 
The solution a(K), A(K), of (4.16) has been tabulated for some values 
of A; further work was required in an extended range for A, for the needs 
of the present paper. (The quantities —a = ¢)(0) and —A = 4,(0) are 
important since they appear in the expressions for the skin friction and the 


rate of heat transfer respectively at 7, = 0.) 


4.3. Flow past the slowly-moving layer of water 
The functions ¢,,%, which are required for the small perturbation 
from ¢y, Yj» satisfy the linear equations 
o, + hy d; o 2¢, ; + dh, >; +o, = (), 
if; +o, dof, +o, 4,9, = 0, ; 
where terms of O(e2) and higher have been neglected. The solution of 
(4.17), chosen so that the total solution satisfies the boundary conditions 


(4.13), is 


1 


(4.17) 


d, ae ld, hy =-a ly’ , 
and the total solution is 


fi = bo—€a7'4,,, 2, = ty —€a ys, (4.18) 


from which the relation between 7* and K is obtained, 


7a 1-K(1- =), (4.19) 


a 


1 
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4.4. Important physical quantities evaluated at the air-water interface 
Before the motion of the water layer is considered it is necessary to 
evaluate the components of velocity, the components of stress, the rate 


ai ar 


of heat transfer, and the temperature, at the air—water interface z = 3 
(or », = 9), since these quantities determine the amount of water produced 
and the manner in which it flows. 

At s=3*, the horizontal and vertical components of velocity are 
respectively 


a” = B, ax, w* = 0; (4.20) 


the horizontal and vertical components of stress are respectively 


. Cw cu) 
. = Baa ae ) 
ox 2 


= xB?2(—a)(p,p,)'241-— 5 (1-K)}, (4.21) 
a’ 
Cc 2 (Cc ow 
v¥ = —p ey pe (+s) 
( Py Ox Cs 
Sp pe 24 By poxt™ + (22, f, + 3h g, ) 0 (4.22) 


(It is not necessary to evaluate the expression in curly brackets explicitly 
unless the distribution of stress along the air—water interface is required. ) 
The rate of transfer of heat is 


and the temperature is 


T* = rf a K(1 <2) (4.24) 
a 


In the above expressions, 8, and 7, are parameters prescribed by the 
potential flow, whereas « and 7* (and therefore, A, a, and A) are quantities 
which must be determined as part of the solution of the whole problem. 

The first shielding effect due to the raising of the air—water interface 
temperature to 7* may be written as the ratio, 


Rate of heat transfer from the air with water layer 





Rate of heat transfer from the air without water layer 
mF uP ig als EP T, 
= K{ —)A(—-]/K(—)4( 
(r-) (7-)/ (7*)4(7*) 
7-7" 2" vi A 
= i Al} / AK? Hl —e— }. 2 
rat Ar.)/ Ar )(-) _— 


When 7, is large compared with 7* and T,,, (4.25) may be approximated 





wn 





m 


by (T,,—T*)/(T,.-—T,,) since A(T/T,) given by (4.16) is very nearly 
constant and «A/a is very small. 

The ratio (7, —7*)/(T,—T,,) is determined from the final solution 
to the whole problem in a later section. 
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5. ‘THE WATER LAYER 
With the expressions (3.4) and (4.20) to (4.24), which determine the 
rate of production and the subsequent motion of the water, we are now in 


a position to consider the equations of motion and the boundary conditions 
for the water layer. 


5.1. The equations of motion 

The equations which govern the motion of, and the transfer of heat in, 
the layer of water of constant density p,, thermal conductivity k,, and 
specific heat c,, but variable viscosity jj, are 


~ 
= 
>I 
ais 
© 
= 
S10 
agin 
> | — 
=) ) 
n> 
a) {— 
4 -~ 
- 
—. 
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YI “~) 
<I 
ee 
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where dissipation has been neglected. 
In order to satisfy the continuity of the expressions (4.20) to (4.24), 
at x = s* we consider a solution of the following form: 
u = Boxf(n2); Ww - (om Be/pe)*fe(ne)s 
r= T,,+(T*-T,,)g2(ne), Nz = (82 pz Mo,)" "3, 
where py, = Hy at T = T,, and 8, is a parameter which must be determined 
from the boundary conditions at the air—water interface s = 3*. A dash 
denotes differentiation with respect to 75. 
With dependent variables f,(y:,) and g.(j,), and the independent 
variable n,, equation (5.1) is satisfied, (5.2) becomes 


and (5.4) becomes 


lhe second momentum equation (5.3) may be integrated to give 


9 CH 


P+ 3 (ps Byx* + py w?) — wy — = constant. (5.8) 
5.2. The boundary conditions 

(i) At the air—water interface 

‘The quantities, u*, w*, 7*, v*, g* and 7* given by the solution of (5.6) 
and (5.7) must be equal to the expressions for these quantities given by 
(4.20) to (4.24); i.e. if the air-water interface is given by 

N2 = 13 = (82/2 Hon)! *3*, 

then at 7, = Nos 


B. f,(n3) = Bre, f2(n*) = 0, (5.9) 
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3 ?(P2 Mam) 


(nf) = (— a)Bi%(a04)*| 1- SK) | (5.10) 


A comparison of terms involving x? in (4.2) and (5.8) gives 


2m 


1p» B3 = Ip., BR. (5.11 
The heat transfer condition is 
317( py fgy,)!2@ £2 mo! (yt) = (— A)B1%(p,4)"® 22K, (5.12) 
Se, te sd | 
and the temperature 7 = 7* gives 
ga(n*) = 1. (5.13) 


The small parameter ¢ is found by elimination of 8/8, from (5.9) and (5.11), 
i.e. 


€ = f(n3)(p.-/p2)*”. (5.14) 


Numerical values of € show that it is O(10~*) so that the neglect of quantities 
of O(c?) in § 4 is justified. 

(ii) At the melting surface 

There are three boundary conditions to be satisfied at the melting 
surface s = 0; these are: 

(a) the mass of water introduced is equal to the amount of the ice lost 

by the body, 

(b) the heat transferred from the water is equal to that transferred to 

the interior of the ice plus the latent heat of melting, 

(c) water is produced at the melting temperature 7’. 

‘The boundary conditions (a) and (4) are most easily formulated in the 
following way. ‘The ice melts at a constant rate w,, so that, for unit surface 
area, in unit time, 

mass of water produced = w,, pg 


and, volume of water produced = w,, p3/p.; 
thus the equivalent velocity of the water at the melting surface is 
t02(0) = ty pa/pa (5.15) 


and since the heat transfer and thus the melting take place in the =-aerection 
and are independent of x, we also have 


u,(0) = 0. (5.16) 


The heat transfer condition is written 


~h( =) = Lw,, ps = Lw.(0)pz., (5.17) 





Oz 


= ¢,(T,,—T_.)@np3 (from (3.4)), 


and since e *) 


we have (=) = [L xe rm i ee a : )]z(0)po. (5.18) 
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Equation (5.18) is a statement that in unit time the heat transfer 
k,(CT ez). 9. causes the mass w,(0)p, (= w,, ps) of ice to be raised through 
a temperature difference T,,— 7, and then supplies the heat to melt it. 


In terms of fy, g. and 7, (5.16) and (5.18) are 





f,(0) = 0, (5.19) 
ts L+¢;(T,, — T 0) ‘ 
g,(0 =m 0», f2(0), 5.20 

g,(0) c(T*—T, Fon fo(0) (5.20) 


and the condition (c) gives 
(0) = 0. (5.21) 


fay 


Equation (5.15) is used to calculate w,, when w,(0) has been determined 
as part ot the solution. 

A solution of the equations (5.6) and (5.7) is required subject to the 
nine boundary conditions (5.9) to (5.13), (5.19) to (5.21); these conditions 
are just sufficient to obtain a unique solution and determine the parameters 
T*, «, n> and f, (for given quantities 7,, 7_,, and 8,) when K, A, and a 
are expressed in terms of 7'* and ¢« through (4.24) and (4.16). 

Clearly, the exact calculations of this system of differential equations 
and boundary conditions would be complicated. In the following section, 
simple approximate relations which give a good overall picture of the 
inter-relation among the various parameters are derived and the results 
are checked by applying a second, more accurate, approximate method. 


6. "THE SHIELDING EFFECT AND THE APPROXIMATE METHODS OF SOLUTION 
‘The important features of the flow in the water layer are the rate of 
convection of mass, which is related to the rate of melting of the ice, and the 
rate of convection of heat which causes a reduction in the rate of heat transfer 
to the ice; this, in turn, controls the rate of melting. In the steady state 
the rate of melting is such that these effects balance. 
The important non-dimensional quantity that determines the rate of 
melting presents itself when (5.20) is written in the form 


f.(0) 





2 - ~—— 9, (0). (6.1) 
~ €3\ is Z 


It is seen that the parameter 
e(T*-T,)  _ p 
oo”: 





is significant in the determination of f,(0) (which is related to the rate of 
melting) in terms of g,(0) (essentially the rate of heat transfer to the ice). 
Although P contains the unknown temperature 7* it is convenient to express 
all quantities in terms of P wherever possible and thus obtain implicit 
relations from which 7* can be eliminated finally. 


A simple expression which describes the shielding effect of the water 
layer is found in the following way. Equation (5.7) is integrated between 








n 
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the limits 0 and 7 to give a heat balance equation 
g3(n*)—95(0) = —o2, | * fe 5 dno 
- 0 
=O, | * fy $2dn2, (6.2) 
U 
since g,(0) = 0 and f,(y*) = 0. Equation (6.1) is written 
, Fam 3 * 7 2 
g,(0) = p> (fo(nz ) —f2(9)). (6.3) 
Thus 
Olas t—SabO) * 
(iz) 82(0) _ 5, p (6.4) 
g,(0) ‘ 
where 
n* 1 ani 
£5 = | fe 82 4y2/\* fe dno (6.5) 
0 /0 


fe. dyy and x | * f; dng are respectively the rates of transfer 


Since x | 
of heat and mass in the x-direction by convection, g, 1s seen to be the mean 
non-dimensional temperature to which the water is raised during convection. 
A measure of the shielding effect is given by the ratio 


Rate of transfer of heat from water to ice 





Rate of transfer of heat from air to water 


2 (0 l 
Sal ’) = = = iK.. (6.6) 
83(y2) 1+%oP i 
The rate of melting may be expressed in terms of g* and 7* by elimination 
of g3(0) from (6.3) and (6.6): 





, 1+g,P 
g.(n*) = —S— %n(—fo(9)), 
I 
which gives 
q* = €3(Tin T f ) t LL T Gy Cal ii 7 Tn JO P3- (6.7) 


Equation (6.7) states that in unit time, for unit area of ice surface, the 
heat g* transferred across the air—water interface first raises the temperature 
of the mass w,, p, of ice from T_, to 7,,,, melts it and raises the temperature 
of the water produced from 77, to the average temperature 7), + g,(7* — T,,) 


during convection. 


6.1. Simple estimates 

Simple expressions may be obtained for the shielding ratio R,, the 
rate of melting, and the thickness of the water layer and the thermal boundary 
layer in the ice. 

We assume that the effect of the motion of the water on the heat transfer 
and stress at the interface are negligible so that quantities of order € are 
ignored when terms of order unity are present. It is also assumed that the 
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main features of the flow do not depend critically on the local velocity and 
temperature profiles which are given the forms 


(2 > (6.5) 


these expressions satisfy the boundary conditions 
f(0)=0, foy*)=0; — go(0) = 0, g,(n*) = 1. 
The average non-dimensional temperature g, is expressed in terms of x, 
by 
and the shielding ratio 
g2(0) I | 

B(n3) 1+%,P 2-2 
so that « is given by the positive root of the equation 

Px? + (12+ 3P)x—12 = 0. (6.9) 
The melting condition (6.1) and the air—water interface stress condition, 
simplified by (6.8), are 





R,, 





~f,:(0) = — = (6.10) 
vr" F2m > 
and 
2 »* B.\32 12 
f,(0) =, 2 = (-a)(3*) (£41) (6.11) 
tn No” Hom By P2 Ham 


which give 


Bo 1/2 : P ,\23 16 ” 13 
= (3) f,{0) = (=) (£41) (- ja!) (6.12) 
M1/ Fam P2 Hom My 
B\!2 .  /Px\13 1/6 - 1/3 ° 
(3) a (=) (fut) (- jaltze | (6.13) 
PD; 7 Fam P2 Kam ; Fe 


Equations (6.12) and (6.13) express f,(0) and 7 in terms of an unknown 
parameter 7* (through P and x). However, the heat transfer condition 
(5.12) is used to express P(7™*) in terms of a basic parameter 
Cp (T.— T,,) 
L+c,(T,,—T_.) 
We have, from (5.12), (4.24), (6.8) and (6.12) 


T -T* ( a 13 (fue 13 ly Hom ) 3 q (= 23 1 4 G5 P 
T*_T, Or, head tees a i 


and since 


and 





S. (6.14) 





T,.-—T* T,-T oe? Pa ee eee 





; a SN: eR: oll “ B+e,7,.-T_.) Pts: 
then S is related to P by 
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The quantity 4 varies only very slowly over a large range of the temperature 
ratio T*/T and an average value is taken. Average values of u3/,,, and a 
are also used. The factor (p,/p2/,,)~'*? depends on the stagnation 
pressure but has been assumed independent of temperature throughout 
this paper. 

Equation (6.15) is used in order to plot quantities of interest as functions 
of S rather than of P. The rate of melting w,, = w(0)p./p, is given in 
non-dimensional form by 


\ 12 P> 23 i 13 .. 13 
tn a ) (=) Pi (Set2e ( _ jalan , (6.16) 
4 By Fam PSN PLes \ Be 





and the thickness of the water layer by 


8 12 Px 13 i 13 % _ \ 2/3 - 
=*(p 1 \ —— (=) ( ta “ Pi (St) = (6.1) / ) 
Hy/ F2m Ms P2\ Pi F1 


When the average values of uF/us,,, A, a, are inserted together with 


om) 


appropriate values of Pir Pa, Pg etc. (for air—water-—ice) we have, 


, 12 
u(t) = 0-96 x 10-7(Px)?3 (6.18) 
by By 
and 
«(Pi B1\"* 7? 
z | = )-22(Pa)}*, (6.19) 
My 


where P is related to S by 
S = 0-24P + 10-91(Px)?9(1 +9, P). (6.20) 
It is seen from (6.18) that the rate of melting w,, is O(10-*) times a 
typical vertical velocity in the air boundary layer, while (6.19) shows that 
the thickness of the water layer is of the same order as that of the air boundary 


layer. 
The thickness @ of the thermal boundary layer in the ice may be defined 
by 
l ° np : k 
aes | (T-T7 jae = th. (6.21) 
Fu Fo Jo W,,, P3 C3 


The non-dimensional form is written 
Qa , \1/2 k Q\12 | 
} J y) 2 lL, Pp 
0( 21) _ ae 1 | - (6.22) 
My P3hifs\ pr J Wy 
When the numerical values are inserted into (6.22) and (6.15) is used, we 





have 
Q 12 
a( =21) = §-69( Pz)-?3. (6.23) 


It is seen that, for ice, the thermal boundary layer is considerably larger 
than the air boundary laver. 

The existence of a thermal boundary layer of this order of magnitude 
show that the present theory is valid for bodies of finite depth during that 


time when the depth is much greater than the boundary layer thickness. 
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[t also explains, in part, why the larger meteors which enter the earth’s 
atmosphere do not become uniformly heated to high temperatures but 
remain relatively cool except in a small region near the melting surface ; 
the thermal boundary layer would also be present if vaporization or burning 
occurred since the character of the solution (3.1) for the temperature 
distribution in the body would be the same. 


6.2. The improved approximation 

As a check on the accuracy of the estimates made in § 6.1 an approximate 
solution of the system of equations (5.6) and (5.7) with the boundary 
conditions (5.9) to (5.13), (5.19) to (5.21) is found by the Pohlhausen method 
for particular values of 7, and T_, anda stagnation pressure of 1 atmosphere. 
In this method e 4 0 so that account is taken of the effect of the motion of 
the water on that of the air. 


Equations (5.6) and (5.7) are written in integrated form 
* 








* - 
[fa gn wad 2 1o 
— fs(n3)- F2(9) +7" = Z | Py dnp. (6.24) 
Hom Ae Ue 0 
9 (n3) —g3(0) = 02, | * fo 82 dno. (6.25) 
The boundary conditions at n. = 75 are 
Q 
P hts B, 
f(ns) = 0 [(ys)=>6 
2 2\2 Bo | 
« 3,\32/ py py \!2 P 
| i i y} . 
—=— f (75 ) a(=) (2 (1 —(1-K)}, . 
Mon, \P2/ P2 Ham a” - (6.26) 
g(y5 ) = 1, 
12 . 
(n*) {3} ES ¢»(T. — T*) 0 
J Y i — ms o~ : 
gee 85 / Po Mom c,(T*-—T,,) o, (A-eA/a) 
where 
12 T* q 
1 _ [ P2 € 
pa(2)°, K=(1--)/(1- S) 
I» p a 
Che boundary conditions at 7, = 0 are 
==) EiOoy= 6 ‘(0) L +e3(T,— T+) #,(0). (6.27) 
r(()) 9 ’ g ee ——aaemner CT > ° (O.2/ 
aad m1 > CI [*- T ) ~ sé 
This system of equations and boundary conditions was solved, with 
T T.. (no conduction into the ice) for a range of values of JT, such 
that 7* took values from 273° K to 373 K, and for the single case 
e eo tee ey he 4500) K by assuming the following profiles 
H2 )3 y 
Is foo +fai | + foo *2 fos ¥3 
No 7] H2 (6 28) 
No nS > 
Lo = S29 + 8a —= + 822-5 + 823 Ge 
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(the number of parameters in (6.28) is one more than the number of relations 
available; the extra condition g"(7>) = 0, from (5.7), was therefore used). 
The results obtained by this method are to be found in $6.3. 


6.3. Discussion of the numerical results 


The first shielding etfect due to displacement of the air by the water 


laver was found in $4.4 as 





Rate of transfer from the air without the water layer | ee i. 
Rate of transfer from the air w ith the water layer we Ai 
approximately. This ratio expressed in terms of S is 
Te-In _1_P_p 
F T* S Cs . 


where P is given in terms of S through (6.20). The parameter P has the 
range 0 < P< 1-25 when the values L = 80 cal/gm, c. = 1 cal/gm °C, 
| adn 100° C and 7), = 0 C are used. ‘The parameter S obtained from 
(6.20) has the range 0 < S < 18-30, the upper limit corresponding to 
T* 100° C and no conduction of heat to the interior of the ice ( i as T ). 
It is seen from the values of Pand S that R, = 1—c,,, c. P Sisalways between 
1-0 and 0-984 so that this shielding effect is a small one (less than 2”,,). 





R Ss Lo 
r Simple Pohlhausen Simple Simple 
estimate 7 T. estimate estimate 
0-125 0-923 ()-923 2-904 0-660 
0-25 QO-859 0-85 4-841 0-654 
0-50 0-756 0-748 8-345 0-644 
0-75 0-675 ()-663 11-710 (0-635 
1-00 ()-612 0-594 15-041 0-627 
1:25 0-568 0-538 18-300 0-620 





Table 1. 


The second shielding effect due to the convection of heat in the layer 


S ° 1 
T Water, given DY 


R 


Rate of transfer of heat from water to ice 





* Rate of transfer of heat from air to water 

s by no means negligible as figure 2 shows. 

The minimum value of R,, corresponding to maximum convection 
f heat by the layer of water, is 0-538 so that approximately 46°, of the 
heat transferred from the air is convected by the water. ‘The good agreement 
between the results of the simple estimate with the values found by the 
Pohlhausen calculation (see table 1) suggests that R, is insensitive to the 
particular velocity and temperature profiles in the water layer; it was found 
that the mean temperature Z, varied little when the profiles were changed. 


F.M. 2L 








522 Leonard Roberts 


The other quantities of interest are the non-dimensional rate of melting 
w,,(p,/u, 8,)'? shown in figure 3, the non-dimensional water layer thickness 
z*(p,8,/4,)'? (figure 4) and the thermal boundary layer thickness 
O(p, 3 fy )§2. (figure 5). 
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Figure 3. Graph of the melting rate. 
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Figure 4. Graph of the water layer thickness. 


‘The rate of melting w,, is O(10~?) times a typical vertical velocity 
(4, 8,/p,)'*? in the air boundary layer, and the thickness of the water laver 
is of the same order as that of the air boundary layer. The good agreement 
between the simple estimates and the Pohlhausen calculation for the melting 


rate w,, 1s due in part to the fact that the quantity (— dayo,, u*)'3, which 
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involves average values, appears in the dominant term of the right-hand side 
of (6.15) and in (6.16) so that when w,,(p,/4, 8,)'? is plotted against S 
the error caused by averaging is unimportant. The agreement between 
the simple estimate and the Pohlhausen calculation for the water layer 
thickness is not so good, however, since the averaged quantity appears as 
(— apy,,'u5)~13 in the expression (6.17) for 2*(p, 8, 4,)'*. 
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Figure 5. Graph of the thickness of the thermal laver in the ice. 
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The simple estimates indicate that the thickness of the water layer 
increases with S through the whole range. However, the Pohlhausen 
calculation, which takes some account of the variation of u#, shows that 


it is possible to have a maximum thickness. Due to the ability of the water 


to convect more easily at the higher temperatures (the viscosity of water 


2L2 
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decreases by a factor of 6-3 in the range 0 C to 100°C), the shielding 
ettect due to convection increases even though the water thickness decreases 
slightly, as S increases. The temperature 7* is found from given values 


of T,, and T_, through 
f* = T,,+Fi1L+¢,7T,,—T_.)i, 


where P is found in terms of 7), and 7T_, through equations (6.20) or 
figure 6. (S is given in terms of 7, and 7'_,.) 
Figure 7 shows typical velocity and temperature profiles in the water 


layer as given by the Pohlhausen calculation. 


nr) 
* 























7. UHE AXISYMMETRIC CASE 
The methods described in the preceding sections may be applied, with 
nly minor amendments, to the axisymmetric stagnation point flows. 


[his is done briefly below. 








The melting of ice in a hot stream of air 525 


7.1. Heat transfer in the ice 
Equations (3.1) to (3.4) remain unaltered since they describe the heat 
transfer in the vertical direction. The solution is 


rf = T_,.+(T,,-—T_. jexp( Saft Ps Cs =), (7.1) 


. 


from which is obtained the rate of heat transfer from the melting surface t 


the interior. 
oT ~ a = 
(A, 7) = W,, p3C3(T,,, — T_,). (7.2) 
g o 


lakes 


7.2. The motion of the air 
At large distances from the body the flow has the potential form 


) > 


u= Bp,Xx, w= —f13, 7 3.5 (7.3 


where x is now the radial distance. 
The non-dimensional variables f,, g, and y, are now given by 





IW 8B. u_\12 . 
aa 3, (71); ~ (2 =) Film); 
P p s 
(7.4) 
E 1 RB. o.\"2 
; +0 P . #E 


The equations of motion and energy for flows with axial symmetry 
become 
SP —-Ati = t8ithi 
£1 T ah g, 0, 
and the boundary conditions are 
0) =0, fKO)=«6 (0) =T*/T., (7.6) 
f, (0) = 1, g,( 2) i 


The system (7.5) and (7.6) with e« = 0 has been solved by the approximate 
method described in $4.2. In the original solution (Levy & Seban 1953) 
(7.5) and (7.6) with « = 0 were shown to represent two-dimensional flow 
past a wedge; these equations together with (7.4) here represent 
axisymmetric stagnation point flow. 

The solution f, = dp, 2; = % has the same form (4.15) but now a and A 


are given by the simultaneous equations 
5K A[2a?— 3(1—K)] = 24a° —5a3[6(1— K)]+ 
+a[Z(1-K)?-3(1-K)], (7.7) 
2445—So,aA? = —}0,(1-K)A. 
The perturbed solution, for small ¢, is again 


fi = dy- s€,, a> by — €a Mb, (7.8) 
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and at the air—water interface we have the following expressions for the 
components of velocity and stress: 


u®* = B, ex, w” = 0, (7.9) 


7t* = xB3/2(—a)(2 Ly oes 1-—K a 
biti bp» Bix? +28, pa{a+ (2g, fit 5h £i)}, m = 0 } 
where 7 satisfies 


w= fi(fi 81) +tMAiga)' + 3 Ay ~ 8 fi oon 


The heat transfer and temperature are given by 


q* = KA( (28, py)" “21 : 


Oo} 


4 > (7.12) 
ry rN €4 j 
re = 7.[ 1-K(1- =) | | 
a } 
7.3. The motion of the water 
The appropriate non-dimensional variables are now given by 
u = Bo xf3(n2); @ = (2ptom Bo/ pe)" *fo(ne); (7.13) 
T= T+ (T*—Ty\ealne) and 19 = (2Bapalttn,)¥*s 


The equations for axisymmetric incompressible flow with variable viscosity 
reduce to 


er sdb pias ) 
——. i ae a 1 = 1 o> he ar _ 
ne ) ‘ ade as Je j (7.14) 


The boundary conditions at the air—water interface are 


fA(nt)=9, Aa)=as 6 
r2 
* | 
~e » feo cys 10 ‘5 ‘ > 
BS *(P2 Ham) * — fo(nz) = (— 4) Bi "(py H1)'? 4 1- a2 (1—K)>, 
Hom <a } 





where 


and at the melting surface 
f,(0) = 90, g.(0) = 0, 
ree: ee ae . > (7.16) 


and 
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7.4. The approximate solutions 
Simple estimates, similar to those obtained for the two-dimensional 
problem, are found by using the profiles 


h = f:0)( 3), 


& =x -= -(l-«) =. 


3 
tw 
a —_~—— — 


The non-dimensional mean temperature g, is unaltered, 1.e. 


ee 
go os = 
52 2 6 


x, 
where 
Px? + (12+3P)x—12 = 0, 

and P has the same meaning as before. The shielding ratio R, is the same 
as that in the two-dimensional problem (equation (6.6)). The rate of 
melting w,,, the water layer thickness, and the thermal boundary layer 
thickness in the ice, are given approximately by 

Wn (Py My B,)' 2 = 1-33 x 10-?(Pa)? , 

3*(p, By/4,)"* = 0-88(Px)*®, 
and 9(p; By/my)'* = 6-26(Pa)** 
(average values of «3, a, and A have again been used). 

The parameter P (which contains 7*) is again related to the basic 
parameter S by 
c 


on py a 13 Piha -13 Hom 13 Pa ** (1435 P). 
Co 2A8 P2 Ham Be Fam ™ 


A second calculation was carried out as a check on the accuracy of the 
results by assuming the velocity and temperature profiles (6.28) and 
applying the Pohlhausen method as described in § 6.2. 


II 


S = 


7.5. Discussion of numerical results 

It is fortuitous that the quantity (a/2A*)'* is very nearly the same as 
the corresponding expression in the two-dimensional problem for most 
of the smaller values of T7*/T,. As T*/T, becomes small and K—1 a 
comparison of (4.16) with (7.7) shows that a/2A% tends to the same limit 
(24/50,)!3 in both cases. 

Consequently the relation between S and P is the same as that for the 
two-dimensional problem to the degree of approximation used here, and 
the relation between S and P remains S = 0-24P+10-91(P«)?*. The 
results for the axisymmetric problem are similar to those for the two- 
dimensional case. The relation between P and S is the same so that the 
shielding ratio R, = (1+, P)~! as a function of S is the same as that for 
the two-dimensional problem. The water layer and thermal layer thick- 
nesses are less and the melting rate is greater than in the corresponding 
two-dimensional flow (the factor involved differs from V2 mainly because 
different average values of a and A were used even though (a/A*)"* has 
the same average values). 
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8. SOME GENERAL CONSIDERATIONS 

‘The flow past a melting body has been considered for the case of a 
semi-infinite body of ice. With this simple geometry the full Navier-Stokes 
equation (with certain restrictions) has been solved in an approximate 
manner. However, if we first make the boundary layer approximation and 
interpret x, z as coordinates measured along and normal to the body we 
again arrive at the ordinary differential equations (4.6) and (4.8); thus 
the results may be applied to any body of large radius of curvature (see, 
for example, Levy & Seban 1953). 

In most applications the high temperature 7’, is produced by a strong 
shock wave and since both y—1 and MW would be small, (WV measured 
behind the shock) terms of O((y—1).M?) would be insignificant. 

Since the rate of melting depends only on the rate of transfer of heat 
and the components of stress at the air—liquid interface (for given conditions 
at 3 + «© and given flow properties), the method may be extended to 
include real gas effects in the air. In this connection Hayes (1956) has 
shown that these effects may be taken into account approximately in such 
a way that the stagnation point flows have similar solutions and the boundary 
layer equations once again reduce to ordinary differential equations. 

‘The more general problem of the melting of a body of finite length has 
no steady solution and the equations of motion reduce to partial differential 
equations; when the body length is large compared with its thermal 
boundary layer thickness the unsteady etfects may be expected to be small 
and some account may be taken of them by a perturbation of the steady 
state solution. 


The author wishes to thank Professor C. C. Lin for suggesting the 
problem, and for his interest and advice during the course of the work. 
The work was sponsored by the Office of Naval Research under contract 
Nonr-—1841(12), NRO62-156 at the Massachusetts Institute of ‘Technology. 
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Energy content and ionization level in an argon 
gas jet heated by a high intensity arc 
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SUMMARY 

A direct current electric arc has been developed as a heating 
device for argon gas. Negligible amounts of electrode material 
are consumed during an operating time of several minutes. Under 
normal operating conditions 50°,, to 90°, of the input electric 
power is transferred directly to the gas. The remaining power is 
absorbed by the water-cooled electrodes. Measurements were 
made to determine the total gas enthalpy and the proportion of the 
enthalpy in directed kinetic energy, random particle motion, and 
ionization energy. From these measurements it is postulated 
that the gas is initially in a non-equilibrium state on leaving the 
arc, but approaches equilibrium relatively quickly when confined 
to a constant diameter jet outside the arc. ‘The gas temperature 
range in these experiments varies from 5000°K to 15 000° K. 


INTRODUCTION 

The study of propulsion systems having a very high specific impulse 
and the investigation of various aerodynamic and aerothermal effects 
occurring at Mach numbers greater than 10 require gas flows with stagnation 
temperatures of several thousand degrees. An electric arc has proved to 
be a very effective means of heating a gas to such temperatures in a 
continuously operating device. Much work has been done in Germany 
in recent years on the water-stabilized are utilizing carbon electrodes. 
The back electrode was usually a rod cooled by the water vortex. The 
front electrode was a carbon plate with an orifice directly in front of the rod 
through which a jet of water and carbon vapour escaped when the arc was 
in operation. The development at the Giannini Research Laboratory of 
a continuously operating arc to heat a jet of inert gas has made it possible 
to study the state of the gas immediately upon leaving the are and at 
subsequent times by confining the gas flow to a water-cooled channel. 
From the gas dynamics of the flow and experimental measurements of 
thrust, mass flow rate, and the total energy balances, it is possible to deduce 
the amount of energy in the gas in the form of directed kinetic energy, 
random particle motion, and ionization. In the results reported here, 
argon was used exclusively. 


* Now at the California Institute of Technology, Pasadena, California. 
8} 
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EXPERIMENTAL EQUIPMENT 

Figure 1 is a schematic diagram of the electric arc used to heat the gas. 
The electrodes consist of two metal plates, one with an orifice through 
which the gas escapes. The plates are water-cooled and enclosed by a 
water-cooled cylinder through which the gas is fed tangentially into the 
arc chamber. The arc strikes between the edge of the orifice and the back 
plate. Gas near the outer edge of the hole receives significantly less energy 
than that which passes through the centre of the arc. The energy exchange 
process appears to be unaffected by either reversal of the polarity of the 
electrode or radial, rather than tangential, introduction of the gas. 
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Figure 1. Schematic drawing of typical device used for arcing through the gas. 
A nozzle is shown in position, electrically insulated from the front electrode 
of the arc 


Figure 2 (plate 1) shows the general arrangement of the equipment 
and essential instrumentation. ‘The thrust developed by the gas jet was 
measured by observing the deflection of a cantilever mounting bar. Flow- 
meters measured the rate of flow of the gas and cooling water. The power 
lost in the electrodes could be determined by noting the temperature 
change of the cooling water. The electrical power input was measured 
during each run and a bank of gauges registered the pressure at various 
points in the arc and confined jet. 


JET CHARACTERISTICS 
Once the gas has left the arc chamber and entered the cylindrical 
channel, diffusion processes tend to equalize enthalpy and velocity in the 
radial direction. Since the Reynolds number is quite large and the nozzle 
length was never more than five times the diameter, boundary-layer effects 
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were ignored in the calculations. Ambipolar diffusion, thermal diffusion, 
and diffusion by concentration gradients are all occurring simultaneously ; 
hence no effort was made to calculate the radial distribution functions. 
However, an estimate of the effect of radial distributions of velocity and 
temperature was made in the following manner. First, Gaussian distribution 
functions were assumed for velocity and temperature. Then the average 
values of these quantities were determined for various distribution parameters 
and these average values compared with the results calculated assuming no 
radial variation. ‘The enthalpy and velocity proved to be least sensitive to 
the radial distribution parameter. With a radial temperature ratio of 20 
(approximately room temperature at the outside) and a velocity ratio 
variation of 4:5, the result of averaging over a Gaussian curve would be to 
increase the velocity and enthalpy by about 20°, over the values obtained 
by assuming no radial variations. ‘The average temperature under these 
conditions would be 70°, higher than the temperature for no radial 
variations. 

A rough estimate of the fraction of the gas affected by the boundary 
layer can be made by measuring the power carried away by the cooling 
water of the nozzle. If r, is the effective radius of the edge of the boundary 
layer, a the radius of the channel, and P,,, the power absorbed in the nozzle 
cooling water, then 


cl 


P,.= 22 | (hy—h)pur dr, (1) 





hy, 
where (4y—/) is the enthalpy decrement across the nozzle. From this 
relation, the mass flow in the boundary layer can be found. This calculation 
assumes a uniform initial distribution of enthalpy throughout the gas, 
which is certainly incorrect, and the answer can be taken as only a rough 
indication of the fraction of the mass flow affected by the boundary layer. 
With these restrictions, the calculations indicate that sometimes as much 
as 60°,, of the mass flow can be within the boundary layer. However, the 
actual physical extent of the boundary layer is a small fraction of the channel 
diameter, due to the much higher gas density in the vicinity of the walls. 
The region of cooled gas can be seen very clearly in figure 3 (plate 1). 

In order to analyse the gas flow, the following assumptions have been 

made: 

(a) electrical energy is transferred to the gas in a channel of the same 
diameter as the hole in the front electrode and mass is added 
continuously along the length of the arc; 

(6) radial variation of temperature, density and velocity are neglected ; 

(c) the energy in electronic excitation is negligible (for the temperature 
and densities used here, Resler, Lin & Kantrowitz (1952) have 
shown that this introduces no appreciable error) ; 

(d) the axial pressure gradient due to the magnetic field of the arc is 
neglected (the magnetic pressure was calculated to be less than 
3°, of the static gas pressure for the current densities and gas 


pressures used). 
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un 


The equation of state for a noble gas in which the temperature is high 
enough to produce ionization may be written 
p= pR( Fa +aT, .s (2) 
and the enthalpy of the gas is given by 
kw oe iw (3) 
Yo—lp 
Definition of the symbols is provided at the end of the paper. The equili- 
brium degree of ionization, x, is given as a function of the temperature 7 
by the Saha equation: 
“* >? AE} A(F®- E}) (4 
a am SS ) 
l—x*p, RI R1 
Values for the free energies (AF°) were obtained from a report by Gilmore 
(1955). By combining these equations with the gas-dynamic equations, 
it is possible to relate the exit gas velocity .,, temperature 7',, and degree 
of ionization x, to the various measured parameters such as thrust F, 
channel exit pressure p,, mass flow rate m, and power input to the gas P. 


\\ hen 


G* 


F>yAp.» 
the flow in the orifice is critical and this criterion was used in the experiments 
to determine whether the flow at the nozzle was sonic or subsonic. 
For the critical flow the following set of equations was derived from 
one-dimensional flow theory (Cann & Ducati 1957). The velocity is 


related to the thrust and mass flow rate by 
te te (5) 
‘Eyl om ; 
The degree of ionization is obtained from the equation 


(Pej) = vk (B+Ab (6) 
ml i) 2(y7-1) m ) 


and the temperature is given by 








= as y F+ Ap, \* e: 
r. Lined ie v0 ; ; 
eee (yo aR ( m ) (/) 
In some cases, the exit pressure was checked by using the result 
F+ Ap, Ap, 
Ap; 4 P ; (5) 


Yot 1 Yor! 

For subcritical flow in the orifice, the above equations are changed to 
the following set. ‘The velocity is now given by 
py, = Fim. (5a) 


The degree of 1onization becomes 


P. _ . FAp a, ae 
a x 70 ( 1+ 70 en | 6s 
( m “i) Yo l m ( 270 =) ( » 
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and the temperature is now obtained from 





sh - FAp ” 
E.-oer). a . 
ett de m*R — 
with the pressure relation given simply by 
Pr = Pa (Sa) 


DISCUSSION OF RESULTS 
(1) ‘Table 1 shows some representative measurements. During the 
experiments the electric input power ranged from 10 kW to 100 kW; 
however, the power in the runs shown in table 1 ranged only between 
30 kW and 50kW. The thrust that the jet developed varied between 
one and ten newtons. 





(1) a . | a 

(2) (3) (5) (6) (/) 
m (4) | : 
h l uz | (ug)q | SI | 
kgm sec : p ip : | 
102 IkW-sec kgm in. 2 m/sec | m/sec | __ sec | 
| | | | 
—_ = ti = “po = Pid mae, en 

2-98 6500 1 0-14 580 | 1800 | 58 
5-07 5130 4 0-12 470 | 1600 | 47 | 
9-06 3390 1 0-08 370} 1300 | 40 | 
| | | | 
5-04 3100 1+. 4-0% 0 1250 | 1240 | 164 | 

ng | 

4-66 4030 L+1-Q* 0 1450 1420 195 





Table 1. Results from a series of tests using argon in the plasma jet. ‘The mass- 
How rate, specific enthalpy, and front electrode length are measured directly. 
The exit flow velocity uz and degree of ionization ag are determined from 
equations (5) or (5a), and (6) or (6a), respectively; while the specific impulse 
is obtained from the definition 5.1. = F'mg; (uz),¢ is the predicted value of 
exit velocity assuming equilibrium in the gas at the exit. A one-inch nozzle 
was added in the two runs marked * in column (3). 


(2) Using the configuration shown in figure 1, without a nozzle, the 
exit velocity obtained from the measured values of thrust and mass flow 
was always much lower than would be expected from a calculation using 
equilibrium thermodynamics. As can be noted from the first three lines 
in table 1, these values differ by a factor of two or three. The ionization 
level of the gas calculated from equation (6) or (6a) was likewise much 
higher than the value predicted for equilibrium ionization, which should be 
negligible for the enthalpy given in column (2). 

(3) The geometrical parameters of electrode spacing, orifice diameter, 
and thickness of the front electrode were varied over a considerable range. 
\lthough some improvement in performance was noted when the width 
of the front electrode was increased, in general the velocity and specific 
impulse remained considerably below the equilibrium values as described 
above. ‘lhe average calculated ionization level remained much higher than 
the expected equilibrium level. 
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(+) The effect of varying the electrical parameters of the arc was also 
studied. The polarity was changed with no marked effect on the gas- 
dynamic performance, despite the fact that the point on the front electrode 
at which the arc struck moved from the outside to the inside when the 
polarity of the front electrode was changed from negative to positive. 
A wide range of average current densities was used with only a second-order 
effect on the exit gas-dynamic properties. There was an indication that at 
higher current densities the deviation from the average equilibrium exit 
conditions was more pronounced. 

(5) The gas enthalpy and pressure were changed over the ranges 
1000 kW sec/kgm to 12000 kW sec/kgm and 18 psi to 150 psi, respectively. 
There was no measurable indication that either of these variables affected 
the noted discrepancy between the calculated and the equilibrium values 
of velocity and ionization level at the exit. 









































T ! | | ii | 
6300 KILOJOULES/KG — 
| | / 
1300 -—+— + ! | 7 + 
1100 : . } 
— 
| [= aoe 
vecocity 99° — T 
AT EXIT IN x | 3300 KILOJOULES/KG 
|e —— 
METERS/SEC_ les 717 | | = 
$00) fo a ® 
S 1700 KILOJOULES/KG 
300] $$ pp ++ 
q j | 
100 l it L l 1 
Vo i 2 
FRONT LENTH OF NOZZLE 
ELECTRODE IN INCHES 


Figure 4. Exit velocity as a function of nozzle length for a 3,16-in, diameter nozzle. 
‘The gas used was argon. 


(6) A series of nozzles was placed over the jet immediately downstream 
from the front electrode, from which they were electrically insulated. 
The nozzles were for the most part cylindrical and of the same diameter 
as the front electrode orifice. The length of the nozzles was varied from 
zero up to 10 nozzle diameters. There was a gradual trend toward 
equilibrium at the exit as the length of the nozzle was increased, and, for 
the longest nozzles, under some operating conditions, the measured exit 
gas velocity and ionization level agreed to within experimental error with 
the calculated equilibrium values. Figure 4 shows the effect of nozzle 
length on the exit velocity of the jet for three ditferent enthalpy levels. 
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It can be seen that the exit velocity reaches approximately a constant value 
for nozzles longer than five nozzle diameters. 

(7) Examination of the back electrode indicated that the arc struck 
over an area equal to that of the orifice in the front electrode; however, 
no estimate could be made of the radial distribution of current density. 
When the back electrode was the anode, the power in watts measured in the 
anode cooling water was numerically about 10 to 15 times the current in 
amperes through the arc. This indicates an anode potential drop of 10 to 
15 volts, which represents the voltage through which the electrons are 
accelerated without suffering collisions with other electrons or gas atoms 
before impinging upon the anode surface. This value agrees quite well 
with values reported by Cobine & Burger (1955) for conventional arcs. 

(8) Attempts to correlate the power absorbed by the front electrode 
and nozzle were never very successful. The best approach appeared to 
be aerodynamic heat-transfer calculations. By using the simple heat-transfer 
equation for turbulent boundary layers, the power absorbed by the cooling 
water of the front electrode could be calculated to within a factor of two. 

(9) The power radiated from the front orifice and the jet was measured 
to determine if any appreciable fraction of the input electrical power was 
lost from the gas in the form of electromagnetic radiation. The spectro- 
grams indicated a fairly intense continuum radiation, as well as numerous 
lines. From calorimetric measurements of the total radiant flux, it was 
found that at most 5°,, of the input power was radiated. This effect was 
not included in the calculations, as it was of the order of the experimental 
errors of measurement. 

The following model has been postulated to explain the behaviour 
noted above. While the gas is passing through the arc, the electric power 
is transferred to a small core of gas at the centre of the orifice, indicated as 
region ab in figure 5; this region, encompassing only a small fraction of the 
total mass flow, does not appreciably increase in radius due to the pinching 
etfect of the arc’s own magnetic field. Once the gas leaves the arc, radial 
diffusion of enthalpy occurs and a larger percentage of the gas is heated up 
(region bc). Simultaneously, for the cases where a nozzle is added after 
the orifice, a cool boundary layer is formed along the outer radius that 
encompasses more and more of the mass flow as the channel length is 
increased. As long as the outward radial flow of enthalpy is faster than the 
radial inward growth of the boundary layer, the flow will behave as indicated 
in (6) above. An axial pressure drop along the nozzle was measured in the 
tests and appears to result from the volume heating effect of the radial 
diffusion of enthalpy and the transformation of ionization energy into random 
particle motion. 

With the above model, equilibrium calculations were carried out to 
see if all experimentally determined quantities could be accounted for. 
When a nozzle was not used, the equilibrium ionization level found from the 
calculations was always low by a factor of two or more from the measured 
value. If the central core of gas in the arc region were in equilibrium, 
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there would be a fairly large radial enthalpy outflow due to the high 
temperature gradients. However, the measurements indicate that no 
appreciable radial enthalpy flow occurs in the presence of the electric field 
since lengthening the front electrode did not affect the thrust for a given 


power and mass-tlow. 
FRONT ELECTRODE 
— ELECTRICAL 


INSULATOR 


_7— NOZZLE 


BACK : 
7— COOL 
ELECTRODE “BOUNDARY LAYER 
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Figure 5. Diagram indicating the proposed energy exchange processes occurring 
in the jet. ‘Che enthalpy is transferred to the gas in the core of the are between 

and 6. From / toc the enthalpy diffuses radially outward until it intersects 

the boundary laver region. From ¢ on there is a hot central core surrounded 


by the expanding cool boundary layer. 


‘Therefore, it is tentatively concluded that the gas is not in temperature 
equilibrium in the are region and in the region immediately downstream 
from the arc; the electron temperature and ionization level in these regions 
1 
i 


must be considerably higher than the equilibrium values. 


The authors would like to extend their thanks to Professor H. W. 
Liepmann of the California Institute of Technology for his many helpful 
suggestions during the course of this work. We would also like to 
acknowledge the invaluable assistance in preparing and proof-reading the 
report rendered by Dr V. H. Blackman of Giannini Research Laboratory. 

This research was supported by the United States Air Force through the 
\ir Force Office of Scientific Research of the Air Research and Development 
Command under Contract No. AF 49(638)-54. 


LIsT OF SYMBOLS 
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a 2. Instrument panel showing the equipment and_ essential instrumentation used during the 
majority of the tests. 


gure 3. Photograph of the plasma jet exhaust showing the region of non-luminous flow near the nozzle 


walls. Standing shock waves can be seen clearly. 
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Energy content and ionization in a hot jet of argon 


} specific ionization energy, R gas constant, 
! front electrode length, e.g. S.I. specific impulse, 
1+1-0 represents a }-in. T temperature, 
thick front electrode with T, electron temperature, 
a 1-0-in. long cylindrical 1, gas temperature, 
nozzle, u gas velocity, 
m gas mass flow rate, % fraction of gas ionized, 
p gas pressure, p gas density, 
p, | atmospheric pressure, Yo 1-67, not necessarily the ratio 
P,, power in the gas, of the specific heats. 


Subscript EF refers to conditions at the channel exit. 


REFERENCES 
Cann, G. L. & Ducati, A. C. 1957 Propulsive properties of high intensity plasma 
jets, Giannini Research Laboratory, Report no. TR-9. 
CosINE, J. D. & BurGcer, E. E. 1955 ¥. Appl. Phys. 26, 895. 


GILMoRE, R. F. 1955 Equilibrium composition and thermodynamic properties of 


air to 24000 K, Rand Corporation, Report no. RM-1543. 
REsLerR, E., Lin, S. & Kantrowi1z, A. 1952 7. Appl. Phys. 23, 1390. 








On the disturbed motion of a plane vortex sheet 
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SUMMARY 


A formal solution to the initial value problem for a plane 
vortex sheet in an inviscid fluid is olftained by transform methods. 
The eigenvalue problem is investigated and the stability criterion 
determined. ‘This criterion is found to be in agreement with that 
obtained previously by Landau (1944), Hatanaka (1949), and 
Pai (1954), all of whom had included spurious eigenvalues in their 
analyses. It is also established that supersonic disturbances may 
be unstable; related investigations in hydrodynamic stability 
have conjectured on this possibility, but the vortex sheet appears 
to afford the first definite example. Finally, an asymptotic 
approximation is developed for the displacement of a vortex sheet 
following a suddenly imposed, spatially periodic velocity. 


1. INTRODUCTION 

The instability of a plane vortex sheet has been investigated for 
incompressible flow in the classical work of Helmholtz, Rayleigh, and 
Kelvin (see Lamb 1945, § 232, § 268) and for compressible flow by Landau 
(1944), Hatanaka (1949)* and Pai (1954). Each of these last three 
investigations ignored the existence of branch points for the eigenvalue 
equation and accepted the eigenvalues given by its two possible branches 
We shall establish the proper treatment of the branch points and the ultimate 
character of the motion by considering an initial value problem for the 
vortex sheet and obtaining an asymptotic solution for large time. Our 
results contirm those of Landau, Hatanaka, and Pai with respect to the 
question of stability but rule out certain of their neutral eigenvalues. 

A second feature of the vortex sheet problem that appears to have 
received insufficient emphasis is the fact that supersonic disturbances, 
i.e. disturbances that have a supersonic wave speed relative to the local 
flow, can be unstable. Such disturbances have often been neglected in 
other stability problems on the (sometimes tacit, sometimes conjectured) 
assumption that they could not be unstable, so that it is especially important 
to ascertain that their presence in the results of Landau, Hatanaka, and Pai 
is not a consequence of their failure to deal properly with the aforementioned 
branch points. We also remark that neutral supersonic disturbances of a 
vortex sheet are significant for acoustic reflection therefrom (Miles 1957). 


* I am indebted to Professor I. Imai for a résumé of Hatanaka’s paper, which is 
not readily available in the United States. 
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A more direct derivation of the eigenvalue equation could be achieved 
by posing a simple, travelling wave disturbance (see (4.11) below). The 
proper restrictions on the eigenvalue equation (and, therefore, the proper 
conclusions regarding the question of stability) could then be inferred 
from Sommerfeld’s finiteness and radiation conditions*, although it is 
to be noted that the propriety of the radiation condition in a stability 
investigation does not appear to have been accepted unequivocally (see 
Lin’s remarks quoted in §6 below, especially his implied question as to 
the nature of a ‘ proper restriction’). Aside from this question, however, 
the possibilities for even asymptotic solutions of initial value problems 
in hydrodynamic stability appear quite limited, and the vortex sheet is 
probably one of the few such problems that is tractable. 


2. FORMULATION OF THE PROBLEM 

We consider (see figure 1) two ideal fluids occupying the half spaces 
y > 0 and y <0, designated by the subscripts + and —, respectively, 
and each characterized by its uniform velocity U _ parallel to the x-axis, 
sonic velocity a_, and density p_. ‘The vortex sheet separating the two 
fluids is subjected to an initial displacement mp(v) and an initial velocity 
n(x); We wish to examine its subsequent motion and, in particular, the 
conditions under which this motion will be bounded (i.e. stable) for large 
time. 








¢ 
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Figure 1. Vortex sheet separating two parallel flows. 








Many of the following equations involve the parameters of only one 
fluid or the other; accordingly, we need include the subscripts + only 
in those equations that involve the parameters of both fluids, with the 
implication that equations devoid of the + subscripts apply to either fluid. 


* I am indebted to Professor G. Carrier for persistently refusing to accept certain 
erroneous conclusions that were improperly inferred frem the radiation condition 
in the original manuscript. Extensive discussions with Professor Carrier then led 
to the attack on the initial value problem. 

+ Stoker (1952, p. 97) has considered an initial value, surface wave problem in order 
to illustrate the validity of the radiation condition in forced oscillation problems, but 
his results do not appear to be directly applicable to eigenvalue problems. 


2M2 
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Assuming small disturbances, the perturbation pressure satisfies the 


wave equation 


D*p 
ay% = ——, (2.1) 
ee Dt? ae 
where D fal U re] >i) 
Dt at Ox ere 


The boundary conditions at the vortex sheet y = n(x, ¢), across which the 
pressure must be continuous and to which the fluid motion must be 


tangential, are 


lo D?n : 
p Piss : cP - a3 y> OL, (Z.3) 
poy Dt i 
while the boundary conditions at infinity are 
lim p < %, lim p < &. (2.4) 


The initial conditions are 


p = 0, t=; y #0 (2-5\) 
Ot ; 
and On ; 
' Ah nD); —— = n(x), C=10, (2.6) 
ct 


We require solutions to (2.1) in y < Oand y > 0 that satisfy equations (2.3) 
to (2.6). 


3. ‘THE TRANSFORM SOLUTION 


We define P, the Fourier—-Laplace transform of p, by 
Ply; m,s)= | edt} e-™*p(x,y,t) dx, (3.1a) 
~ 0 - x 


1 f -ytia : 
“—t a abies dm | e P(y; m,s) ds, (3.1 b) 

47771 s . ins ‘ 

where y is a positive real number such that all singularities of P lie in 

A\s\ <y. Similarly, we define N(m,s) as the Fourier—Laplace transform 





P(x, ¥5 1) 


ix 


of n(x, t) and 


Ny(m) = | e-™*ng(x) de, (3.2a) 


Vo(m) = | | é mf ig) + us) dx, (3.2 b) 


dx 
as the Fourier transforms of the initial displacement and _ transverse 
velocities imparted to the fluids. We also introduce the operator 
S =s+imU, (3.3 
corresponding to the substantial time derivative of (2.2). The transform 
of the wave equation (2.1), subject to the initial conditions (2.5), then 


reads 
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while the transforms of the boundary conditions (2.3), subject to the initial 
conditions (2.6), read 


| oe = P ” _ = = —= S2N _ SNo- Ve (3.5) 


Solving (3.4) subject to (3.5) at y = 0, we obtain 


imp, p_(U,—U_)(S,V,_+S_Vo.)e*=" 











? ge . -_ 3.¢ 
P_(y; m,s) are eT (3.6) 
and N(m, s) = (p. by )(S., No +I 0 ) + (p_/m \(S No+ I 0 ) (3.7) 
0,+0 
where 
O = pS?/u (3.8) 
and 
pe = {m* + (S/a)*)?, Abu > 0 (3.9) 
js} 
4 1{s} 
—(a—U)(ma— im) # R{s} 4 i(Imla — mu) 
“= Uns ta 


-1(Imla+ mU) 


Figure 2. Cuts for p(s, m) in the s-plane for Figure 3. Cuts for p(s, 7) in 
U < a,m, > 0) me > 90. the s-plane for m real. 
The condition .#\u} > 0, which is a consequence of the second 
finiteness condition (2.4), will be satisfied everywhere in the complex s- and 
m-planes if we choose the branch cuts from the branch points 


stim(U+a)=0 (3.10) 
to be lines on which -#/! = 0 or, equivalently, nu? < 0; the latter condition 
yields 

a’m, m, + (s;— Um ,)(s,+ Um,) = 0, (3.11a) 
a*(m= — m5) + (s, — Um,)? — (s, + Um,)? < 0, (3.11 b) 
where 


m = m,+imMs, S = $,+ 15>. (3.12) 








542 John W. Miles 


Equations (3.1la, b) define the cuts as segments of a hyperbola in the 
m-plane for fixed s, or in the s-plane for fixed m; if m, = 0 these segments 
degenerate into straight lines in the s-plane, and similarly in the m-plane 
if s,= 0. The cuts in the s-plane for U <a, m, > 0, and m, > 0 are 
illustrated in figure 2; if m, <0 the cuts from (U+a)(m,—im,) and 
—(a—U)(m,—1im,) must be asymptotic to Um,+i1x and Um,—ix. The 
cuts for m, = 0, as in the problem to be considered in the following section, 
are shown in figure 3. 

To be sure, the condition A{u.} > 0 need not be satisfied everywhere 
in the s- and m-planes but only along the paths of integration for the inverse 
transforms of P. when +y > 0. We imply here only that it is sufficient 
to choose the cuts defined by (3.11a,b); in the final analyses they may 
be deformed in any manner that ensures convergence of the integrals 
and the satisfaction of the second finiteness condition (2.4) by p.. 


4. SPATIALLY PERIODIC DISTURBANCES 
It suffices, for the question of stability, to consider a spatially periodic 
(fixed m) disturbance. Such a disturbance will result if we assume both mp 
and m, to be periodic in x, but the algebra is simplest and the results equally 
illuminating for the special case 
n(x) = 0, tte(x) = v,¢'*7, (4.1) 
where vy is the amplitude of an initial, transverse velocity imparted to the 
vortex sheet at ¢ = 0 and « the wave-number (figure 1). We must take x 
to be real in order to satisfy (2.4), and we may assume it to be positive without 
loss of generality. The required transforms of (4.1) are 
N, = 0, V = 277, 8(m— x), (4.2) 
where 5(m— x) denotes the Dirac delta function. Substituting (4.2) in (3.6) 
and (3.7), writing the inverse transforms according to (3.1 b) and carrying 
out the m-integrations yields 














p40. OU ey OO 8, +8 et ee ; 

. ne it = 4.3 
P+(% st) 2a | ws ps p(Q,+Q_) ia 

and n(x,t) = x | A | & 0 ae E ‘| esttiar do. (4.4) 


where m = x in p, S, and O. 
It is found convenient, especially in relating our results to more 
conventional, travelling wave analyses, to introduce the change of variable 


s= —12C, y = ae, (4.5) 


where c is a complex wave speed. Writing also 





; ee FA 1/2 
8 -£=[(' =) -1] , SB > 0, (4.6) 
4 a 


and F = 1Q/x = p(c— U)?/B, (4.7) 
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(4.3) and (4.4) become 
p.p(U_—U_)zvy s*>"€ (2c—U,—U_)eis(4-4+8 21) de 











pi(x,y,t) = | 
. : = 
2a | ee 8B, B_(F.+F_) 
(4.8) 
nd 
to (27 [(p./B.)+(o_/B 
n(x,t) = — = 2) (p 3) eilz—t) de. (4.9) 
LUE } etic F_+F.- 
lhe branch points and cuts for 8 are shown in figure 4, as also are its phase 
angles on the axes #'c! = 0 and Aic! = U (the latter is significant for an 


observer moving with the fluid, i.e. in the coordinates x— Ut and y). 


I{s} 


| |(w/2) 


(7/2) (O) R{S} 


PITIFIFIFIII AS 


(0) -(0-U) | (U a+U (7) 


~— - 


(7) | 


FIFIFIIIIIIIIIIIT 





' 
| 
Figure 4. Cuts for 8 in c-plane for m real; the phases of 8 on the real axis and on 
# ic} = U are indicated in parentheses. 


There are a total of four branch points, namely, c = U.+a_., U,-—a,, 
for the integrands of (4.8) and (4.9). We may assume U_ > U_ without 
loss of generality (if U_ > U_ we have only to replace y by — y), and there 
are then four possible configurations of the branch points, corresponding 
to (see figure 5) 

*, 0< U.-—U_<a_-a. 


\ 
A,: 0< U_-U_<a_-a_ 
\ , 


The three possibilities A,,, are not significantly different (in particular, 
as shown in § 5, the essential character of the eigenvalues is the same for all), 
but in B we find that it is impossible for a disturbance to be subsonic, 
i.e. Ajc}—U <a, with respect to both fluids. We also note that the 
cuts for 8. may be deformed freely in evaluating p., as given by (4.8); 
moreover, only two finite segments of the real axis appear as cuts for the 
integrand of (4.9), the phase jump of the quantity in square brackets 
vanishing if 8. and §_ are either both real or both imaginary. 
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We now consider the nature of the elementary, exponential disturbances 
from which (4.8) is synthesized. ‘This may be done conveniently in the 
coordinates x— Ut and y, especially if we introduce the transformation 

c= U+asec6, fp = tan§, (4.10) 


under which 





exp{iz(x—ct+P y\)} = exp{ixsec 6[(x— Ut)cos#+ ysiné—at]}. (4.11) 
If @ is real, as it is for the points on cuts in the c-plane, this corresponds 
to an outgoing (incoming) plane wave of sound in x—Ut and y for 

















R{C} 
-(a_-U_) (04+U,) 
-(a,-U,) (a_+U_) 
A\0<U,.—U_<e_-—<6, 
R{c} 
—(a4- U4) (a_+U_) 
~(e..~U.) (a4+U4) 
Ao: O<U4—U_<04—90- 
R{C} 
—(o.—U_} (a_+U_) 
— (04-U4) (a4+U4) 
Ax la,—a_I<U,—-U_<a,+a_ 
R{c} 





—(a--U_) (U4—a4) 
(a-+U_) (U,+a+) 
B Ua SUS 6,5F 6 


Figure 5. Cuts for 8. and f_ in c-plane for m real. 

sinf >(<)0; more generally, points in -7{c} > (<)0, correspond to 
outgoing (incoming) waves that fall off exponentially as y — a. Of 
course, the total disturbance is confined to y < at, for if y > at we may 
close the path of integration in .7{c} > € to obtain 

p(x, y, t) =0, v >at. (4.12) 
It follows that the envelope of the total disturbance is the outgoing wave- 
front y =aat, as is otherwise directly evident from the definition of a. 


5. ‘THE EIGENVALUE PROBLEM 
The vortex sheet at y = 0 is stable only if the integrand of (4.9) has 
no poles in .Z{c} > 0; otherwise it is unstable. These poles correspond 
to the zeros of 


F .(c)+F_(c) = p_(¢ —~U_)?/B.4 p_(c— U_)?/B (5.1 a) 


=p a (B + B-1) + p a- (8_+ B=") (5.1 b) 





wn 
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and may be identified as the eigenvalues for elementary disturbances of 
the type (4.11) subject to the boundary conditions (2.3) and (2.4). We 
remark that (5.1) does not contain the wave-number and is therefore 
significant for any disturbance that has the phase velocity c. 

We may simplify the algebra considerably, without unduly restricting 
the results, by assuming equal specific heat ratios in the two fluids. Then, 
in virtue of the requirement of equal static pressures, 


p.a* = p_a*, (5.2) 
and we may factor (5.1 b) to obtain* 
(8, + B_)(8.B_+1)=0. (5.3 a) 
Alternatively, the trigonometric substitution (4.10) yields 
sin(@, + @_)cos(#@, —@_) = 0. (5.3 b) 
We consider first the zeros of 8, +8. Remarking that these zeros 


must be real in consequence of the finiteness condition .7{f} > 0, and 
that 8. and B_ can be of opposite sign only on the real axis interval 
U_+a_<ec< U,.-—a, (which arises only for U,—U_>a,+a_), we 
deduce that 


ae ate = (8), yy ~t opie, Of 


I {c} 


a hy 


*\ 1{8,2.} 








(a) (b) 

Figure 6. Mapping of the upper half c-plane on the 8 ,8 _-plane. 
This zero corresponds to 6.+6 = 7 in (5.3b). It clearly represents a 
neutral, supersonic disturbance and was obtained by Landau and Pai 
along with the spurious, real zero of 8, —B_ = 0. 

We may investigate the possible zeros of 8 8 _+1 (equivalent to those 
of cos(#. —@_)) in-4Z}c} > 0 by applying Cauchy’s principle of the argument 

* The eigenvalue equation considered by Landau and Pai was 
(8, B2. —1)(82 —B2) = 0, 


while Hatanaka considered (82 82 —1) = 0. 
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to the contour consisting of the real axis of the c-plane, indented above the 
branch points of 8. and 8_, anda semi-circle of radius tending to infinity 
in the upper half-plane. ‘This contour and its map on a 8,8 -plane (the 
Cauchy~Nyquist diagram) are shown in figures 6 (a), (6) on the assumption 
{ U’_ >a_+a_; the general shape of the 8,8 ~ contour is similar for 
the other branch point configurations of figure 5, only the order in which 
the points are traversed changing. It is evident that the 8,8 contour 
will encircle the point — 1 once if (8.8), > —1 (as shown in figure 6 (6)) 
or will pass through —1 twice if (8.8_), < —1, from which we infer 
that 8.8 +1 must have either one zero in .Z{c! > 0 or two zeros on the 
real c-axis, respectively. Similarly, we infer from a contour indented 
below the branch points and closed in the lower half of the c-plane that 
8.8 +1 may have one zero in.Z{c} < Oif (8. B_),) > —1. We also remark 
that both of the 8.8 contours pass through the point + 1 twice (indepen- 
dently of the location of the point (8.8_),), corresponding to the two 


real zeros of 8.8 —1; these cannot be equal, being less than U_—a 
and greater than U_ +a_., respectively*. 
The critical point (8.8 _)) = —1 may be determined by requiring 


8.8 +1 and its derivative with respect to ¢ to vanish simultaneously ; 
now, however, having proved that 8,8 —1 cannot have a double root, 
we may work with £2 8% —1 or, more conveniently, 


"- (- ) ' (+) at (5.5) 


which is readily shown to have the double root 


oa (= gE SB ) (5.6) 





a*? + a’ 
at 
(U,, — U_)}*® = a48 + a2, (5.6b) 
We conclude, therefore, that 
(8,.B_),»2 —-1 as (U,—U_)<S (a2 +a2%)3?, (5.7) 
The limiting case U, = U_ (= U, say) requires special consideration, 


since the two complex conjugate zeros of 8.8 +1 then coalesce at ¢ = U. 
We infer that an infinitesimal, tangential discontinuity in a compressible 
fluid will exhibit dimear (rather than exponential) instability, thereby 
generalizing Rayleigh’s result (Lamb 1945, § 232) for incompressible flow. 
The foregoing results are summarized in table 1, and we conclude 
that, subject to the restriction p,a* = p_a*,, a necessary and sufficient 
condition for the stability of an inviscid vortex sheet with respect to small 

disturbances is 
Cae i > (a? 34 gi/3)3 2 (5.8) 


* These last zeros were included in the analyses of Landau, Hatanaka, and Pai; 
indeed, the chief flaw in these analyses, insofar as they aim only at a stability criterion, 
is their failure to prove that the complex zeros of £2 B? 


B,B_+1 alone. 


—1 may be charged to 
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in agreement with the predictions of Landau (1944), Hatanaka (1949), 
and Pai (1954)*. 





Zeros of (5.3) 























Relative speed range 
P Fici > 0 
ic} = O ic 0 
| * (unstable) aii a 
a —_—__—___— | — ——— 
O U_-l | 2* 0 0 
——EE ae a —_ — — —_— } | ——__— a 
)} < |U,—U_| < a.+@ | 0 l | 
| 
a ee ee —— a ee 
a,ta_<|U.,-l[ | < (a + q?!8)3/2 | 1 1 1 
(a2? +a28)3?2 < |U,—U_| 3 0 0 
| | 








Table 1. The possible eigenvalues of equation (5.3). 


* double root atc = U. 
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Figure 7. The stability boundaries of (5.12) and (5.13). 


The explicit determination of the zeros of 8.8 +1 requires the solution 
of a quartic equation, but in the special case of equal sonic velocities 
(a. = a_=a) this may be reduced to a quadratic equation having the 
roots 

c= 1(U,+U_)+ fa[M?+4-4(M?4+1)2]}'2, M=(U,-U_)/a, (5.9) 
* Neither Hatanaka nor Pai gave the result (5.8) in explicit form, but their final 
results are in agreement therewith. 
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where M is the Mach number of relative flow, which, according to (5.7), 
must exceed 2%? for stability. 

A special case of greater practical significance than that of equal sonic 
velocities is that of equal stagnation enthalpies (e.g. the vortex sheets in 
shock interaction patterns) where 








LEIS icky we GEE ewe, (5.10) 
. yl : y-1 

Introducing the notation 
M,, = U.,/a,, m=a.j/a., (5.11) 


and combining (5.10) with (5.6b), we obtain the parametric representation 
M_ = (y—1)71(1 + m?3)-32(1 — m2) — $(1 + m?3)82, (5.12 a) 
M, = (y—1)-'m—1(1 + m?3)-32(1 — m2) + 3m—1(1 + m?3)3 2, (5.12 b) 


This stability boundary is plotted in figure 7 for y = 1-4, as also is 


that obtained by combining (5.10) with U, —-U_ = a,+a_ (see Lin (1953) 
and § 6 below), namely 
M_ = 3[(3-y)—(v+1)m]/(y-1), (5.13a 
M, = $[-—(3-y)+(y+1)m“]/(y-1). (5.13 b) 


6. THE ROLE OF SUPERSONIC DISTURBANCES 

Supersonic disturbances, defined as those for which Afc}—U >a, 
have played a somewhat ambiguous role in previous investigations of 
hydrodynamic stability*. Referring to these disturbances for a boundary 
layer, Lin (1955, pp. 70, 71) states that: 

‘‘One would then expect the disturbance outside the boundary layer 
to be a wave with non-diminishing amplitude at infinity instead of an 
exponential decay of the amplitude. ‘There is no discrete characteristic 
value problem for such disturbances unless some proper restriction is 
imposed. In fact, these ‘supersonic disturbances’ have not yet been 
fully studied. 

In all the existing theoretical analyses of the stability of the boundary 
layer in a gas, supersonic disturbances are assumed to be insignificant. 
This is based on the conjecture that the energy associated with such 
disturbances would propagate from the boundary layer in the nature of 
acoustic waves. Additional theoretical work and experimental evidence 
on this point are highly desirable.”’ 

If we had presumed the impossibility of unstable supersonic distur- 
bances for our model of a vortex sheet, the condition 


jU,—-U_| >a,+a (6.1) 
*If | pic} -U| > a the elementary disturbance of (4.11) has a supersonic phase 
velocity, Aic} — U, along the coordinate axis x— Ut. The phase velocity with respect 


to the fluid is [1+(#{8})*] -1"(Aic} — U); if c and f are both real, the latter velocity 
is equal to the sonic velocity a, as is directly evident from (4.11). 
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would have appeared a priori as sufficient for stability. The more severe 
criterion (5.8) is a consequence of the existence of unstable supersonic 
disturbances for 
a.+a_<|U,-U_| < (a3 +a?3)8?, 

This, in our opinion, implies the possibility of their existence in related 
problems of hydrodynamic stability and adds emphasis to the concluding 
sentence in the above quotation from Lin. In particular, we should expect 
(5.8), which constitutes a necessary condition for the stability of a mixing 
region of zero thickness, to be significant for the stability of a laminar 
mixing region of finite thickness with respect to inviscid disturbances 
having wavelengths large compared with the thickness. The latter problem 
has been studied by Lin, who suggested the criterion (6.1) as sufficient 
for stability on the tentative hypothesis that only subsonic disturbances 
could be unstable. It is not clear that the analysis presented here is directly 
comparable with Lin’s*, but the present results certainly enhance the 
desirability of its extension to include supersonic disturbances. 

We remark that the presence of a boundary at some finite value of y, 
such as the wall in the boundary layer problem or the virtual boundaries 
for symmetric and antisymmetric disturbances of a symmetric jet (cf. Pai 
1951), may be decisive for the stability of supersonic disturbances. 
A vortex sheet parallel to a boundary at which either p, = 0 or p = 0 is 
considered in the Appendix and is found not to admit neutral, supersonic 
disturbances; it does not necessarily follow that unstable, supersonic 
disturbances are impossible for this configuration, but it seems likely that 
this also should be so. 


7. ASYMPTOTIC EVALUATION OF ” FOR @_ = @, 


We now return to the solution obtained in §4. An explicit evaluation 
of the integrals (4.8) and (4.9) in terms of known functions does not appear 
to be feasible, but an asymptotic development for n is straightforward. 
We may further simplify this development, without losing any significant 
features, by assuming 

p. =p. = p, a,=a_=4, Gi) 


so that (4.9) reduces to 
Uo -O+ie e! x( r —Ct) dc 


rn ee es 


*In applying boundary conditions, Lin assumes that the displacement of any 
streamline must be small compared with the thickness of the mixing region, so that 
the approximations appear rather different than those adopted here. It is possible 
that the approximation of small displacement could be relaxed to one of small stream- 
line slope simply by an implicit change of variable, e.g. a von Mises transformation 
in which y is replaced by the stream function. A second feature of the finite mixing 
region that is absent in the model of a vortex sheet is the inner viscous region (Lin 1955, 
p. 136), where c = U and the inviscid differential equation has a singularity. The 
role of this singularity as the thickness of the mixing region tends to zero would 
require special attention. 





n(x,t) = — (7.2) 


2raa? J 
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The symmetry of this special case may be emphasized by shifting to a 
coordinate system (x ,y) that moves with the mean velocity cy of the two 
fluids, where 











Ne = X—Cpf, @= HU, +0U_). (7.3) 
We also introduce the dimensionless wave speed w, such that : 
C= Cy+ au, B.. = {(wz 4M)? —1}2, c=“, (7.4) | 
where M is defined by (5.9); (7.2) then becomes 
Up etrto perie e iaatw dw —: | 
n(x,t) = — 2 | eee eae 73) 
Ge es ( 


The poles of the integrand are given by (5.9) as 

w= +i, A = [(M?+1)!?- (4M?+4+1)]??. (7.6) 
The complex w-plane, the original path of integration C, and the cuts and 
poles of 8.8 +1 are shown in figure 8; if W > 23? the poles lie on the 
real axis between + }M—1, rather than on the imaginary axis. 
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Figure 8. The z-plane for (7.5). 


If ¢ > 0 we may deform C into the four, clockwise loops C, 53, shown 
in figure 8, plus a contour at infinity in the lower half-plane that makes a 
null contribution to the result. The phases of §.8_ on the top and bottom 
of C, are }z and $7 and conversely for C,; hence, the integrals over the 
top of C, and the bottom of C, are found to differ only in the sign of the 
exponent, and similarly for those over the bottom of C; and the top of C,. 
Utilizing these data to combine the four integrals over the tops and bottoms 
of the cuts and evaluating the contributions of the poles with the aid of 
Cauchy’s residue theorem yields 
vpe'*" { sinh(xAat) 

wa \A(M?+1)"? 


2. ' M+ sin(xatw) dw 


a(x,t) = 








Suwa 1 ife® = (EM — 1PM 1 
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If M > 232, 
sinh(xAat)/A = sin(« Ajat)/ A}. (7.8) 
The imaginary part of the integrand of (7.7) vanishes like w—|}M+1) 1? 
at the end points and has no other singularities or points of stationary phase. 
Accordingly, integrating once by parts and making use of 2-8 (11) in 
Erdélyi’s (1956) monograph, we obtain the asymptotic approximation 
n(x,t) = ve" ( sinh(xAat) _ 
wa \AUM? +1)? * 


? ay 1/2 
+ ae) (xat)~3 *eos| (4. + 1)xat F 





I 


a 


|+ O[(zat)-*2]\, (7.9) 
where the sum of the two terms corresponding to the upper and lower 
choices of sign is to be taken. 

We remark that a similar development is possible for the pressures, 
but then (assuming the cuts as in figure 5) the integral will extend along 
the entire real axis, and there will be contributions from points of stationary 
phase; moreover, if 7 > 2 the contributions of the poles on the real axis 
(one at w = Oif 2 < M < 23? or three at w = 0 and +/A if M > 23?) must 
first be separated out. 

APPENDIX 
Vortex sheet near boundary 

We remarked in § 6 that the presence of a boundary at some finite value 
of y may be of decisive importance with respect to the stability of supersonic 
disturbances. We attempt to throw further light on this point for the 





vortex sheet problem by placing a wall at y = —h, so that 

Ua Of. -h<y<®B, (Ala) 
- U,, y>0. (Alb) 

We then have, in addition to the boundary conditions (2.3), 
ph, = 6, y= —-h, (A2) 
Assuming the elementary disturbance of (4.11) and imposing (2.3) 
and (A2), the resulting eigenvalue equation is found to be (cf. (5.1a)) 
p < y+? _ en ae ee (A3) 
If the boundary condition p, = 0 is replaced by p = 0 at y = —A, it is 


necessary only to replace coth by tanh in (A3). 

We remark that: (i) B=! coth(—7iz8_ A) is a single-valued function of c, 
so that the only branch points of (A3) are atc = U,+a,; (it) (A3) reduces 
to (5.1a) as h ow unless f_ is real (supersonic, neutral disturbance in 
y <0), so that (A2) differs from a radiation condition as h-> o only if 
B_ is real; (iii) there can be no supersonic, neutral disturbances in the 
upper medium, since f. is real for such disturbances, while the second 
term in (A3) is imaginary for all real c; (iv) the statements (i) to (iii) 
remain valid if the hyperbolic cotangent is replaced by the hyperbolic 
tangent. 
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We surmise from the foregoing remarks that supersonic, neutral 
disturbances may not exist for a boundary layer near a fixed wall or for a 
symmetric jet (for either symmetric or antisymmetric disturbances). 
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